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A method for numerical integration 
on an automatic computer 
By 
C. W. CLENSHAW and A. R. CURTIS 


Abstract 


A new method for the numerical integration of a ‘‘well-behaved” function over 
a finite range of argument is described. It consists essentially of expanding the 
integrand in a series of Chebyshev polynomials, and integrating this series term 
by term. Illustrative examples are given, and the method is compared with the 
most commonly-used alternatives, namely Simpson’s rule and the method of Gauss. 


1. Introduction 


We consider in this paper the basic problem in numerical integration, the 
evaluation of the integral of a function /(x) known numerically within a finite 
range of integration aSx<b. For the user of an automatic computer, it is 
desirable that a program to perform this task should give the result to a specified 
accuracy in a wide range of circumstances, and give an indication of the reduced 
accuracy achieved when these circumstances do not obtain. 


Finite-difference methods are, in general, poorly suited to automatic com- 
putation because their storage requirement is large. Lagrangian formulae, which 
involve only function values, are usually preferred. The two most commonly 
used are probably Simpson’s rule and Gauss’ formula. The former has the 
attraction of simplicity of form, with convenient binary weights, while the latter 
minimises the number of ordinates needed. 


We propose here a method based on the term-by-term integration of the 
expansion of /(x) in Chebyshev polynomials. A unique advantage of the 
method is that its accuracy may be checked easily before the integration is 
completed. It also shares some of'the advantages of the other two methods 
mentioned; as with Simpson’s rule, the number of the ordinates may be doubled 
if necessary without previous work being wasted, while considerable economy 
in the number of ordinates is also achieved. Moreover, the method will provide, 


with little extra complication, values of the indefinite integral f f(x) dx through- 


a 
out the range (a,b). In contrast, Simpson’s rule yields values of the indefinite 
integral at tabular points only, while Gauss’ formula is quite unsuited to in- 
definite integration. 

Although we are primarily concerned with the integration of a non-singular 
function in a finite range, we may nevertheless observe that an infinite range 
may frequently be transformed to a finite range, or approximated by a large 
finite range. Further, some integrands with weak singularities are amenable 
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to treatment with the proposed method. In general, however, special types of 
integral such as 


fe-* f(x) dx and Te-* f(x) dx, 
0 —oo 


are probably best evaluated by special methods such as the Laguerre-Gauss 
and Hermite-Gauss formulae respectively (see, for example, HILDEBRAND [1] or 
KopPaAL [2]). The second of these integrals can also be evaluated efficiently by 
using the simple summation formula of GOODWIN [3]. For integrals of oscillating 
functions over an infinite range, LONGMAN [4] has described a simple method 
based on the Euler transformation of series. 


2. The methods of SIMPSON and GAUSS 


Srmpson’s rule is given by 


Jt(x) dx GhH(0) + 4f(h) +f(2A)}. (1) 


This is the three-point ‘‘Newton-Cotes” formula. The corresponding N-point 
formula is obtained by integrating the interpolation polynomial of degree N — 1 
which is equal to the integrand /(x) at each of N equally-spaced values of the 
argument x. In practice, we subdivide a given range (a, 6) into an even number 
2m, say, of equal segments, and use the summed form of (1), namely 


6 
fre dx (b — @) (fo + 4h + 2fe+4fat+ 2fatees +4 fom—1t+ fem)» (2) 


where 





2m 


i= t{ (2m—r) are), 


Similar summation procedures may be adopted with the other Newton-Cotes 
formulae. In general, the formulae involving an even number of segments are 
preferred, since each of them has the same order of accuracy as the formula 
involving one more segment. Of these preferred formulae, SIMpson’s rule is the 
simplest, with weights which are powers of two. The higher order formulae 
have coefficients which must be stored, and moreover these coefficients are large 
and vary in sign when N is large, thereby magnifying the efiect of rounding 
errors. 

If (2) is to yield a satisfactory result, m must be made large enough to ensure 
that the truncation error is negligible. This may imply that many function 
values have to be evaluated; for functions which can be rapidly computed, 
however, this is usually considered worthwhile in order to permit the use of a 
simple and flexible formula. In such cases, however, in order to deal with functions 
having regions of rapidly changing behaviour, a programmer often finds it con- 
venient to allow the interval to vary as the integration proceeds from step to 
step. The variation is usually controlled by comparing at a given step the ap- 
proximate integrals obtained using two intervals of integration, # and 2h, say. 
If these results agree, then an interval 4h is tried, and so on until disagreement 
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is found, whereupon the last agreed value is accepted, and the corresponding 
intervals are used to start the next step. If the first two results disagree, the 
interval is repeatedly halved until agreement is reached. 

Comparison of two results is commonly used to assess the accuracy of an 
integral evaluated on an automatic computer by Simpson’s rule applied at a 
constant interval. However, this check is unreliable, as may be shown by a 
simple example. If 


f(x) =%%coshx—cosx, a=—1, b=4, 


the values obtained from (2) are 0.47955 46 for m=1, and 0.47955 51 for m= 2, 
whereas the true value of the integral is 0.47942 82. 

The possibility of check failure can be greatly reduced by comparing at least 
three successive results; although the extent of the computation is thereby in- 
creased, this seems a highly desirable safeguard for a general quadrature program. 


The Gauss integration formula (HILDEBRAND [1] or Kopat [2]) is given by 


6 n 
S Hx) dx (b—a) YH, f(x), (3) 


where 

x, = 3(d — @) t,+$(b+a), 

1 

iil (1 —@7) {Fy (t,)}? ’ 
and ¢,, ¢,,...,¢, are the zeros of the Legendre polynomial P,(¢). Formula (3) is 
exact if f(x) is a polynomial of degree not exceeding 2" —1. This considerable 
power is achieved at the expense of introducing a set of irrational arguments 
and weights, whose values need to be stored. 

The checking of the Gauss formula is more difficult than that of SImpson’s 
rule. Repetition with a new m requires a completely new set of abscissae x, 
and weights H,, and we cannot make use of previously computed ordinates if 
n is increased. The most commonly used check for an isolated integration is to 
compare the results of applying the Gauss formula at and at +1 points. That 
this is not a completely satisfactory procedure is shown by another example; 
the Gauss three- and four-point formulae give 1.585026 and 1.585060 respectively 
for the integral 





+1 
dx 
whereas its correct value is 1.582233. 

However, there is an important class of problems for which the use of the 
Gauss formula is advantageous. if it is necessary to integrate over a given 
range several functions whose behaviour within the range is similar, it may be 
possible to find an m adequate for all cases by a thorough investigation of a 
small sample. 


3. Method of Chebyshev expansion 


The proposed method derives its power from the economy with which a function 
can be represented by its expansion in Chebyshev polynomials, and the ease 
15* 
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with which such a series can be integrated. If {(x) is continuous and of bounded 
variation in (a, b), then it can be expanded in the form 


f(x) =F) =$4+a4T,)+aR()+-:-, (aSxSbd), (4) 
where 


T, (t) =cos(rcos4#), t= are ; (5) 


(see N. B. S. [5]). On integration we have 


<7) fioax=froa=an+6.n0 +b, T,(f) +>, (6) 


where 


b, = SS (7 = 1,2...) (7) 


(see CLENSHAW [6)). 
The value of b) is determined by the lower limit of integration; thus 


by = 20, — 20, + 20,—-:-. (8) 
The definite integral is given by 


h +1 
pia LIMP ATE Th tht 
=2(,+0,+6;+---), (9) 


while the indefinite integral is given by the sum of the series (6), the evaluation 
of which is described later. 


The coefficients in the expansion (4) may be calculated after first observing 
that any polynomial of degree N in x may be written in the form 


f(x) rnp epty enc eacene cs inne iterss ia besten ton (10) 


-r"6 T(t), (—1<#S1). 


Here >'”’ denotes a finite sum whose first and last terms are to be halved. 
The coefficients in (10) are then given by 


2 z nmuvs 
=H 2, F, cos 7 ae (11) 
where 
F= =F (cos 47). 


This is a consequence of the orthogonality of the cosine function with respect 
to the points ¢,= cos ie , expressed by the equations 
N 0 (¢+7) 
"cos =*S cos 7IS =] N (t=j7=0 or N)}. (12) 


N N 
ia 4N (¢=7+0 or N) 
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Finally, any function which satisfies the conditions necessary for its Cheby- 
shev expansion to converge can be approximated to any required accuracy by 
a finite series of the form (10) the coefficients of which are given by (11). 


It is important to observe that (11) can be written in the form 
2 e uy 
=H 2 F, T, (¢,), t, = COS >. (13) 
This form is similar to that of (10), so that the same machine subroutine can be 
used to evaluate the right-hand sides of these relations. Such a subroutine may 


be based on the recurrence method described by CLENSHAW [7], which we now 
outline briefly. 


The series 
F(t) = $4) +4, T,(t) +--- +4, T, (d) (14) 
is evaluated by calculating successively the numbers ¢,,,c,_1,...,¢ 9 from the 
relation 
C, = 26,41 —C,494+4,, (7Sn), (15) 
starting with c,,.;=C,,2,=0. The sum of the series is then 
F(t) = (Co — ¢2), (16) 


This result may be verified by substituting for the a, in (14) their expressions 
in terms of the c, from (15) and using the recurrence relation for the Chebyshev 
polynomials 

T,41(¢) — 2t7,( + T_1(t) =0. (17) 


It only remains to consider the determination of the value of N. We propose 
that if no information about the behaviour of the integrand is available, the 
Chebyshev coefficients be evaluated using an arbitrary initial value of N in 
the formula (11). We then examine the coefficients of the Chebyshev poly- 
nomials of highest order in the integrated series; if these are small we may 
assume the chosen N to be adequate. We consider the smallness of three con- 
secutive coefficients of a convergent series to be a reliable indication of the 
negligibility of al/ higher order coefficients. For instance, in the case of indefinite 
integration we might require by, kby_, and k?by_, to be negligible, where & is 
an arbitrary constant in the range (0,1). The value of k could be taken to be 
unity; however, a smaller value will often save considerable computing time 
while only slightly increasing the chance of check failure. We suggest the use 
of k=}. 

In the evaluation of a definite integral alternate terms vanish, and we suggest 
as a criterion that three successive non-zero coefficients should be small. This 
“‘safe’’ number of three may, of course, be increased if we desire to reduce even 
further the possibility, already remote, of an erroneous result. 

If by, Rby_, and k?by_. are not all negligible, we replace N by 2N and again 
compute the a, and b, of high order. The previously computed values of F, 
will again be used; they are the alternate members of the new set 


F, = F (cos aw) (s=0,1,...,2N). (18) 
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It is convenient to take N=4 first, since no smaller value can give an acceptable 
number of negligible coefficients, and then to proceed to N=8, N=16 and so 
on, to an upper limit N=M, dictated by the requirements of the programmer 
and the characteristics of the computer. 


4. The case of slow convergence 


The only situation remaining for consideration is that in which the desired 
accuracy is not attained at N=M. In such a case, we seek to determine the 
reduced accuracy achieved. 

We observe that if the integrand may be expanded in an infinite Cheby- 


shev series 
co 


F(t) = 349+ Ai Tt) + 4B) +: = 2 AT, (19) 


then the coefficients a, given by (11) are related to the A, by 





2e ars < mis 
ur ‘ ¢ 
a=5 p> cos NW 2 A, cos Ny 


ni A, > Aey-, a Ay, + Aqn-, > Agyiy 2 Agy-- ed (20) 


41 
The approximate value of f F(t) dt, calculated from the coefficients a, is 
-1 
2a, 2% i, 2an—2 an (21) 


fy = &— St 3.5 (N—3)(N—1)  (N—1)(N+1)’ 





with N even, whereas the correct value of this integral is 





I me hg ~ 243... NE dine : (22) 


Assuming that the coefficients A, may be neglected when r=>3N, we may ap- 
proximate the error Iy —I by 








Ey = 2Agy— Fy (Aaw-s + Aeyy2) —-*°— WOH (Ay+2+Asy—2) + 
= 2An+ r 
+ pa (N+2r—1) (N-F2r+1) ° (23) 


Since the series (19) converges for —1S¢S1, the coefficients A, satisfy an 
inequality of the form 
|4,/s**, (2M), (24) 


where Ky is independent of 7, and is assumed to be as small as (24) will allow. 
Then we have 


1 1 
|Aen-es + Awn+as| Ss Ky (yy -f wes) 


N K 4st 
= Ky a < Sy (1+ Sy) (s<}NQ). 


(25) 
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It is clear that the coefficient of Agy;2, in (23) does not exceed ———__?____ 
in absolute value for s=1, 2,..., } N—1. Therefore (2s—1) (2s+1) 


Eel <Bu fy +) ty (t+ aa) 





(26) 


4N-1 
pate x fhe 1 2 4s? 2 


The largest of 2|ay_,|, 2|ay_2| and |ay| may conveniently be used as an 
approximation to 2Ky/N. Thus, if the upper limit N=M is reached before 
the desired accuracy is attained, we may assess the magnitude of the error to 
be that of the largest of these quantities, or, more conveniently and almost 
equivalently, the largest of 4N |by43+y-1+6y_3|, 4N |by4it+by_,| and 
4N |by41| (see (7)). 

In the case of an indefinite integral similar considerations apply, but here 
all the coefficients 6, are required, and the estimate of 2Ky/N may therefore 
be based on the values of by,,, by and by_, instead of by,,, by_, and by_s. 
This procedure is illustrated in the second of the examples which follow. 





5. Examples 


As our first example we consider again the definite integral, 
+1 
° dx 

#44 2240.9 ’ 


and suppose we require the result to six-decimal accuracy. Taking N=4, we 
evaluate the integrand at x,=cos} 2s (s=0O, 1, 2,3, 4) and, applying formula 
(11) with r=4, obtain a, and hence b,=+0.00609544, which is not negligible 
to six decimal places. We therefore take N=8, evaluate the integrand at the 
intermediate points and apply (11) with y=8 to obtain a,, and hence b,= 
— 0.00007 329, which is still not negligible. Again doubling N we find a new 
@,, and b,,=-+0.00000002 which can be neglected. We therefore calculate a,,, 
keeping N= 16, and thence 20,,, and, finding this to be only +0.00000004 we 
proceed to a, and gz6,3;=—0.00000007. This also is sufficiently small, so we 
accept N = 16 as being large enough, and proceed to compute the terms 3, (7 odd), 
shown in Table 1. 


Finally, the integral is estimated by the present method to be 
2 (by + by + 0; + +++) = 1.58223 296, 


which is correct to eight decimals. 
As an example of indefinite integration which also illustrates the problem 
of slow convergence we consider 


S\x +4hlidx. 


The coefficients b, produced with N=4, 8 and 16 are unacceptably large, and a 
computer program would continue to double N until sufficient accuracy was 
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obtained, or until it reached its limit M. For a program working to about 30 
binary figures, M=64 might be chosen, but for the purpose of illustration we 
take M=16. We find for the coefficients 6, in the integrated series the values 
given in Table 2 


























Table 1 Table 2 
’ by r b, r by 
1 + 0.85844 113 1 +0.707670 9 + 0.000062 
3 — 0.07354558 2 +0.127592 10 — 0.001338 
5 + 0.00645 162 3 + 0.020533 11 +0.001061 
7 —0.00015279 4 — 0.022044 12 — 0.000180 
9 “—0.00010230 5 + 0.008786 13 — 0.000427 
11 + 0.00002 844 6 +0.001172 14 + 0.000516 
13 —0.00000436 7 — 0.004192 15 — 0.000161 
15 + 0.00000030 8 + 0.002548 16 — 0.0001 78 
17 + 0.00000 002 17 +0.000118 


From (8) we obtain b)= 1.250724, and from Table 2 we find that the values of 
64 | big + by¢+ 445], 64 |b,7+,6| and 64 |b,,| are 0.014144, 0.003840 and 0.007552 
respectively. Hence we conclude that the indefinite integra] may contain errors 
up to 0.0141. As a check, the approximate value of the definite integral is 


+1 
S\x+4lidx=4b,+6,+---=1.466900, 
-1 


compared with the correct value 1.460447 ,and is thus in error by about 0.0065. 


It may be of interest to note that a program based on the proposed method 
has been prepared for the DEUCE at the National Physical Laboratory. The 
upper ae M of the number of terms is 64, and the program evaluated the 


integral J |x+3lbax with an error of 0.00078. This may be compared with 


the error ~ 0.00317 committed by the 32-point Gauss formula, and of 0.00036 
by the 64-point Gauss formula. We see that the Chebyshev formula, which 
is much more convenient than the Gauss, may sometimes nevettheless be of 
comparable accuracy. 


Conclusion 


The principal advantages of the Simpson and Gauss methods are to som 
extent shared with the present method of Chebyshev expansion. 

Srmpson’s rule has simple and convenient weights, and its accuracy is easy 
to check because the interval may be halved without wasting the previous 
computations. Its main disadvantage is that many ordinates may be required. 


The method of Gauss requires the least possible number of ordinates. Its 
most serious drawback is that this economy may be lost in checking the accuracy 
of a result, because the weights and abscissae of the m-point formulae are com- 
pletely different for each value of mn. 
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The present method effects a compromise. It is as easy to check as SIMPSON’S 
rule, yet achieves considerable economy in the number of ordinates required, 
and will readily provide values of an indefinite, as well as a definite, integral. 
This combination of economy and reliability make it a very suitable basis for 
general subroutines for numerical integration on an automatic computer. 


The work described above has been carried out as part of the research programme 
of the National Physical Laboratory and is published by permission of the Director 
of the Laboratory. 
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Uber die asymptotisch beste Approximation stetiger 
Funktionen mit Hilfe von Bernstein-Polynomen 
Von 


C. G. ESSEEN 


1. Einleitung 


In der Approximationstheorie fiir reelle Funktionen haben die sog. Bernstein- 
Polynome eine gewisse Anwendung gefunden, indem man z.B. den WeierstraB- 
schen Approximationssatz mit ihrer Hilfe sehr leicht beweisen kann. Ist f(x) 
eine im abgeschlossenen Intervall [0,1] stetige Funktion, so ist das /(x) zuge- 
ordnete Bernstein-Polynom B,,(x) n-ter Ordnung durch 


B(x) = °(2) Pal), (1) 
definiert, wobei 
Pno(x) = (") 271 — a)". 
Es sei w(6) der Stetighkeitsmodul von f(x), d.h 


w(8) = max |f(x)—f(9)| 05051 x,yE [0,4]. (2) 


Es ist bekannt, daB es eine Konstante K gibt, so daB 
\f(x) — B,(x)| SKo(n-*) n=1,2,.... (3) 


Es sei fiir ein festes 





%, = sup max \f(¥) —B, (*)| 
f(x) OSx51 © w(n~4) 


und x= sup x,. Kiirzlich hat SIkKEMA [1] bewiesen, daB 
n=1,2,... 


1% < 1,093 785 (4) 


und dadurch friihere Ergebnisse von Popoviciu [2] und LorENTz [3; S. 20] ver- 
scharft. SIKKEMA zeigte auch, daB 


lim max It(#)—Bn()| <i+ ai = 1,053991. (5) 


n—>oo 0S*51 w(n~*) 


In dieser Arbeit werden wir den asymptotischen Approximationsgrad (5) naher 
studieren. Das Hauptergebnis ist die Bestimmung einer Konstante C, so daB 


lim max lf (¥)— Bn (4)| <<. (6) 





n—>co OSx51 w(n~4) 
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Dabei ist C in dem Sinne optimal, daB man jedem ¢>0 ein n=n(e) und eine 
Funktion /,(x) mit dem Stetigkeitsmodul w, (6) zuordnen kann, so daB 


max |/,(x) — B,(x)| 2 (C — e)@,(n-'). 


Os*sl 
Ist 
x 2 
O(x) = + [ e *adt 
()= 7 fe (7) 
die Normalverteilung, so ist C durch 
C=2))(v +1) (O(2v + 2) — D(2r)) = 1,045 564 (8) 
»=0 
bestimmt. Es gilt also bs 
C= lim x,. (9) 


Bei der Bestimmung von C werden wir auch einige andere Ergebnisse hinsichtlich 
des asymptotischen Approximationsgrades erhalten. 


2. Einige Hilfssatze 


Wir erinnern zuerst an einige Eigenschaften des Stetigkeitsmoduls. Ist w (6) 
durch (2) definiert, so ist w (6) fiir OS 6S eine stetige nicht abnehmende Funktion 
und w(0)=0. Ferner ist w(A6)S(A+1) w(6) fir OS A6S1. Man sieht auch 
leicht ein, daB der Stetigkeitsmodul von w (6) identisch mit der Funktion @ (6) ist. 

Um die folgende Darstellung nicht unterbrechen zu miissen, geben wir in 
diesem Abschnitt einige Hilfssatze an. 


Hilfssatz 1. Ist p,,,(x) durch (1) defintert, so gilt 


aaa" oo (10) 
wy | Bu -an 


wobet k eine beliebige positive Konstante ist. 


Die Ungleichung (10) ist nichts anderes als die Tschebyscheffsche Ungleichung 
der Wahrscheinlichkeitsrechnung, angewandt auf die Binomialverteilung. Siehe 
auch LORENTZ [3; S. 6]. 


Hilfssatz 2. Es sei OS 22 \n x(4 —x). Dann gilt 


p> Puy (x) S2e™. (41) 


Iz al vos |/ =4 
" n 


Diese Ungleichung wurde von S. BERNSTEIN [4] bewiesen. Siehe LORENTZ 
[3; S. 18]. 

Unter ¢(m) verstehen wir von nun an eine nicht naher bestimmte Funktion, 
die durch eine Funktion K(m) majoriert werden kann, 


|e(n)| SK(n), n=1,2,.... 


» 
—-—x 
n 
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Dabei ist K(n) eine beschrankte, von x, /(«) und w (6) unabhangige Funktion mit 


lim K(n) =0. 
Hilfssatz 3. Fiir —2 lognsts<2 log n, n-*<x%<1—n-§ und y=nx+ 
t |/n x(1— x) 1st (v—nx)? 
Puy (%) ‘7 oN owl. 2 e eeti-e . (1 2) 


2a x(1—x) 
Bekanntlich erhalt man die Relation (12), wenn man die Approximation der 
Binomialverteilung durch die Normalverteilung betrachtet. Hierbei ist x in der 
Regel konstant und ¢ in einem beschrankten Intervall variabel. Wir wollen 
aber die Relation auch in den Fallen anwenden, in denen sowohl ¢ als auch x 
zwischen den Schranken des obigen Satzes variieren. Es ist jedoch nicht schwer 
zu beweisen, daB die Relation (12) auch dann giiltig ist, z.B. mit "Tilfe der aus- 
fiihrlichen Darstellung in FRECHET [5; S. 103]. 


3. Abschatzung des asymptotischen Approximationsgrades 


In diesem und dem folgenden Abschnitt sei f(x) eine auf [0, 1] stetige Funktion 
mit dem Stetigkeitsmodul w (6). Wir setzen 


A, (x) =|f(*) — B,(x)|. (13) 
Dann ist P e 
v=0 v=0 











Es hat sich als zweckmaBig erwiesen, die Falle, in denen x sehr nahe 0 oder 1 
liegt, gesondert zu betrachten. 


Satz 1. Ist OS x<n-* oder 1—n-*'<x<1, so gilt 





A,, (x) < w(n-*) + 2n-* w(n—*). (15) 
Aus (14) erhalten wir 
A,,(x) = dol io *|) Pn» (2) = 2; + 2p. 
v=0 . 


° ° - . v 
In 2, summieren wir iiber alle » mit | -—% 
n 


<ni jt” , in 2, iiber alle » 
>ni ata. Nun ist 





. v 
mit |— — x 
n 








2150 (ns 285") < w(nin-* n-) = w(n-*), 


Aus Hilfssatz 1 ergibt sich ferner, daB 
2, Sw (1) n-FS (m2 +1) w(n-2) n-F < 2n-* ew (n-*%). 
Damit ist (15) bewiesen. 


Satz 2. a) Es set 
n-§oxs1—n-, (16) 
Dann gilt 


2V/log n : 1? 
A, (x) < v5 [ ole Yer 2dt+|e(n)|w(n-4). (17) 
—2)log n 
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b) Es seten w(d) ein gegebener Stetigkeitsmodul, n>4 und x eine Konstante 
am Intervall (16). Dann gibt es eine Funktion f(y, x) des Argumentes y mit dem 
Stetigkeitsmodul w (6), OS 6S4, so daB umgekehrt 


2Vlog n 
f(x, x) — B,(x)] = im f ola /eajer —|e(n)}o(n-3). (48) 
ee 


Beweis. Aus (14) erhalten wir 











Ay(*) S So (|2 — x|)yo(2) =2i + Za, (19) 
wobei me 
| 2/2) log in 5, (20) 








Iz - x| >2 | ee in 2. 

n n 

Wir fangen mit der Abschatzung von 2, an. Aus Hilfssatz 2 folgt fiir z= log n 
2, S0(1) 2 <2 V+ a (nt) = e(n)w(n-d, (21) 

Die Bedingung z= |/log nx \x(1—x)n ist wegen (16) sicher erfiillt. 


In 2; ersetzen wir p,,(x) durch die Abschatzung (12) des Hilfssatzes 3 und 
erhalten 


“(5 --) 
(v—n x)* 
415 Mi) Dia ———"— g 8mx(1—z) | 


Eine elementare Rechnung zeigt, daB wegen (16) und (20) 


- 2-2. a ee 
eg 2nx(l—z) S(1+|e(n)|) fe @nz(l—z) dy 


und daher gilt 
v+1 
ols) fate 
a1 (1 e Sas0—2) dy. 22 
+ |e(n)| l) Qa eae J2xnx(1—x) J an 
‘ etl ‘ _ __(w—nx)* 
In (22) ersetzen wir jedes Glied durch Sy er w (| co x|} e %nx(1—2) dy, 


Der so begangene Fehler ist héchstens gleich ” 


Vea =f Jo(lz — o(|% — |) 


n (23) ist zwar 
u 
o(=-*)-#(s-*) 


(u—n x)* 


e 2mx(1-2) dy. (23) 








so(>), (24) 
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aber wir kénnen diese Ungleichung nicht ohne weiteres anwenden, da w (5) von 
n 


derselben GréBenordnung wie o (7) sein kann /z.B. wenn w(6)= wnnrel } 





n log — 
Gilt fiir » die Ungleichung (20) so variiert o (| _ x|) zwischen 0 und 





@ pena eee ) S (2 |/x(4— x) logn+2) w(n-!). Die Anzahl der Inter- 


valle in Z, der Lange ' , fiir die die Zunahme von w(6) groBer als n-* w (n-4) 
ist, ist also héchstens 
N = 2n* (2 //x(41 — x) logn + 2). 


Diese Intervalle geben dann einen Beitrag zu dem Fehler (23), der nach (24) 
nicht gréBer ist als 





1 N vil is 
Der Beitrag der iibrigen Intervalle ist kleiner als 


(un—x)? 


—Hey(n-4) fe 88-2) dy = - o 
” “ asenaca | egal aot ) 


Zusammenfassend ist also (23) kleiner als | e(#)| @(n—4). SchitieBlich erhalten wir 
nach einer Variabeltransformation 


2 Viog » beading 
2, 5 (1+| €(n)|) Vi | ole yj) e 2dt+|e(n)|w(m-4). (25) 
—2 log» 
Nun ist 
2 log» saiactaanans atin p 
1 *(1—z)\_-FZ 
yaa f ole ye 2 dt 
—2 log» 
1 a \t| . 
“tax J lay)? *# 
—2)logn 
2 4 2 Viog » a 
, ned Is d 
spi([e fre T Gta 
< 2 (w(n-4) (@(2) — G0) + w(2n-4) (G(4) — H(2)) + ---) 
< w(n-4) 2 (G(2) — G0) + 2(H(4) — H(2)) + ---). 
Also folgt 
2Vlogn ——————— e 
1 *(1—%*)\ - oni 
as pw y= )e 2 dt<Co(n-4), (26) 
~2Viogn 


wobei (x) durch (7) und C durch (8) definiert sind. Aus (19), (24), (25) und 
(26) folgt die Richtigkeit des ersten Teiles des Satzes. 
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Um die Umkehrung dieses Ergebnisses zu beweisen, betrachten wir die Funk- 


tion 
f(y, x) =@(|y— x]), 


wobei x konstant ist und im Intervalle (16) liegt. Man sieht unmittelbar ein, 
daB der Stetigkeitsmodul von /(y, x) wenigstens fiir O<d6<4 mit w(6) identisch 
ist. Betrachtet man den Beweis des ersten Teiles des Satzes genauer, so sieht 
man leicht ein, daB 
2 Vlog» rn) 
f(x, x) — B,(x)| => Yas [ w (le 24) e *dt—|e(n)|w(n-l), nB4. 
—2Vlogn 
Damit ist dann auch der zweite Teil des Satzes bewiesen. 
Aus Satz 1 und 2 erhalten wir den folgenden 


Zusatz. Ist w(n—*)=0(w(n-*)), n> 0, so gilt 


2 Viog n 
wage) Bec ae= fol aye)e t+ olor. 
—2Vlogn 


Beispiel. Es sei w(6)=K 6*, wobei K und « Konstanten sind und 0<«<1. 
Dann ist 





: att 
lim max Baba Ser ee I il e © 2 a 
nc 05451 ey (n~4) Viz 2 Va 


Aus dem zweiten Teile des Satzes 2 folgt, daB diese Konstante nicht verbessert 
werden kann, wenn man alle Funktionen /(x) betrachtet, die zur Klasse Lip (a) 
geh6ren. 

Bemerkung. Es scheint notwendig zu sein, die Falle gesondert zu betrachten, 
in denen x sehr nahe 0 oder 1 liegt. Ist z.B. x= -, so wird die Binomialverteilung 


nicht langer durch die Normalverteilung, sondern durch die Poisson-Verteilung 
approximiert und der Grad der Annadherung mit Bernstein-Polynomen ist 


<= Konst. w (5) : 
n 
4. Bestimmung von lim *» 
nm—> Oo 


In diesem Abschnitt werden wir den folgenden Satz beweisen. 
Satz 3. 


lim max /f()—Bl)l <¢, 
n—oo 05x51 w(n —) 
wobet 
C=2)> (y+ 1) (O(2v + 2) — D(2r)) = 1,045 564 
y=0 
und 
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Die Konstante C kann nicht verbessert werden, da man jedem e>0 ein n=n(e) 
und eine Funktion f,,(x) mit dem Stetigkeitsmodul w,,(d) zuordnen kann derart, daB 


max \fn(*)— By (*)| ye all 
0sx51 Wy, (n~ Tie 


Beweis. Der erste Teil des Satzes folgt unmittelbar aus den Satzen 1 und 2 
und aus (26). 


Um den zweiten Teil zu beweisen, definieren wir /,(x) folgendermaBen: 


0 x= 2 
1 1 1 1. ; - 
f,(x) = 44 = te <*s * + ya a eine positive Konstante 





, 1 1 1 
linear 3 2*s7+T 
In dem Rest des Intervalls [%, 1] ist /, 4 rekursiv definiert und zwar durch 
h(¥+ 7) =a+tf(); +4+2 sest4 Att, 
V Ir n 
bee 4:2. HY 
In [0, $] ist f(x) dadurch bestimmt, daB die Funktion symmetrisch beziiglich 


x= ist. Der Stetigkeitsmodul w, (6) von /,,(x) ist offensichtlich identisch gleich 
/,(6) fir OS dS} und 








w,(n—*)=a, nZQ4. (27) 


Ohne Einschrankung kénnen wir annehmen, daB m eine gerade Zahl ist, 
n=2m. Dann ist 


Ban ($) — fam (5) = 3; Peme(5) ome ahs) =2_De Pome ($) om (a) 


v=0 v=m+1 
m+ [2m] m+ [2/2 m] , 
. 2a 2s bom (a )+ a uae i me (2) * sil 
m+ (3/2m] ‘ 
’ sh a Penr(3) + 2 “> 


wobei R ein Fehlerglied ist, das von den Punkten — stammt, die in den Inter- 


k 1 


vallen 0< on _ liegen. (Bemerkung: Eine nahere Unter- 














m \2m 4(2m)?* 
suchung zeigt, daB fiir ein irrationales n gilt 
v k Bil. -_ , mes \n 
— Te > a5: y=0,1,2,...; =1,2,. ve F 


Es kommen also tatsachlich keine solchen Punkte vor.) Aus (11) und (12) erhalten 
wir folgende Abschatzung fiir R: 
Konst 


Rs 24D) Pam taltalimss(3) a“ (29) 





Uber die asymptotisch beste Approximation 213 


Wenn wir nun die £3,,,(}) in (28) durch (12) des Hilfssatzes 3 approximieren, 
erhalten wir ahnlich wie beim Beweise des Satzes 2: 


Ba m() — fam(¥) = 20 





Eo + 1) (G(2» + 2) — G(2y))} — e(2m)| —R_ (30) 


und aus (30), (29) und (27): 
|Bam(E) — fom (4)| = (C — e(2m)) w((2m)-4). 


Damit ist der zweite Teil des Satzes bewiesen. 
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Optimal Points for Numerical Differentiation * 
By 
HERBERT E. SALZER ** 


Abstract 


New n-point rth derivative Lagrangian numerical differentiation formulas employ 
the best irregular locations of points, A. From the standpoint of highest degree 
accuracy for derivatives at a single fixed point (mth degree accuracy proven to be 
the highest exactly attainable for any rv). B. From the Tschebyscheff standpoint of 
minimal largest |remainder| over an argument range. In B the dominant term in 
the remainder is minimal for arguments at the zeros of r*® order integrals of Tscheby- 
scheff polynomials specialized by addition of suitable (ry —1)t® degree polynomials 
chosen to produce real, distinct locations of points within or fairly close to the range 
of optimization. First and second derivative formulas up to nine-point, are obtained 
with remainder estimates. 


Introduction 
It is fairly well-known in both interpolation and numerical integration that 
there are certain advantages in accuracy in choosing arguments that are not 
regularly spaced, but at special irregularly spaced — Thus in #-point La- 
grangian interpolation choice of the _— x%;,t=1,...,, at the zeros of the 
Tschebyscheff polynomials 7,*(x) = sQr— 8) 2, 





; cos {m arc cos x}, i.e., at x;=Cos 

minimizes the largest absolute vane “of the factor I (x —x;) in the remainder 
n t=1 

TT] (x —%,) # (0)/n! over the range (—1,1) in x***. Similarly, for evaluating a 

i=1 


definite integral by numerical integration, we have Gaussian-type quadrature 
formulas where the irregularly located points are at the zeros of certain orthogonal 
polynomials (e.g., LEGENDRE, LAGUERRE, HERMITE, depending upon the range 
and weight function of the integral). The resulting -point formulas are exact 
for any (2% — 1)" degree polynomial in the integrand. Furthermore by integration 
of the 2-point Hermite osculatory interpolation formula in terms of /(x;) and 
f'(x,;) at say, x; equal to the zeros of the Legendre polynomials, for finite limits 
of integration, the coefficients of /’(x;) vanish (by orthogonality) and the factor 





[{fte-%9 } dx intheremainder /”(@) {LT (<—x)}*as 
1 


1 “sl (2m)! 4 Gel 


is minimal. 





* Presented at the Eleventh International Congress of Mathematicians Edinburgh, 
Scotland, August 14—21, 1958. 
** Staff Specialist, Flight Performance and Analysis Group, Convair (Astro- 
nautics) Div. of General Dynamics Corp., San Diego, 12. Calif. 
*** We assume a linear transformation of variable so that the new * ranges over 
(— 1, 1) when the original x ranges over (a, b). This will be made more precise later on. 
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The purpose of this present note is to investigate the existence of optimal 
locations of points in numerical differentiation formulas which would give greater 
accuracy than the usually employed equally spaced arguments. The term ‘“‘best”’ 
is employed in the following two senses: 


A. The first sense refers to numerical differentiation at a fixed point only 
in which we seek to find m-point differentiation formulas that express a particular 
order derivative as a linear combination of the function at certain points, which 
shall be exact for polynomials of degree higher than » — 1, the degree of accuracy 
attained in the usual formulas obtained by differentiation of the n-point Lagran- 
gian interpolation polynomial. This will be analogous to the Gaussian-type 
quadrature formulas for definite integrals where n-point formulas attain (2m — 1)" 
degree accuracy. It turns out here that Gaussian-type numerical differentiation 
formulas can give some improvement in degree, but considerably less than in 
corresponding quadrature formulas. 


B. The second sense, which will be dealt with more extensively here, refers 
to that of minimizing the largest absolute value of the remainder when the 
Lagrangian numerical differentiation formula is applied for a variable argument 
lying within some definite range. In other words, the problem is to find the 
most favorable irregular location of the m fixed arguments x; in the formula 


= (n) 
(1) he Oy Hed + Rule), 
where L(x) is the Lagrange coefficient [J] (x —x;) / IT (x%;—~%;), for some 
jt; 7=1 jt; 7=1 

range of x that will reduce |R,,(x)| in some Tschebyscheff sense for x within 
some interval (a, 6). This problem has a satisfactory answer (where the choice 
of the points x; is independent of the /(x)) by considering only the dominant 
part of R,(x), where appreciable improvement is possible. 

This present study of optimal irregularly located points for numerical differ- 
entiation is believed to be new and not related to some published remarks of 
other writers, all of which have been concerned only with the size of the interval 
on the basis of just equally spaced arguments (e.g. [7], pp. 104—110). 

Applications of these results are extremely numerous and widespread, due 
to the fundamental importance of the basic operation of numerical differentiation, 
and will not be cited here except for the following use in astronautic and missile 
science: Suppose one were observing positions of a moving object at different 
times and wished to compute from them either the velocity or acceleration, with 
minimal truncating error. The results here would specify the optimum times at 
which to make the observations. Similarly there might be observed functions 
of position (e.g. tilt or pitch) whose rates of increase with respect to distance 
might be needed, and again the new formulas would specify the most favorable 
locations of the observation stations. 


A. Optimal Points for Differentiation ata Fixed Point 
We may propose the problem of finding distinct points x;, +=1, 2,..., ”, 
such that in some formula resembling (1), (which for amy location of x; is always 
16* 
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exact when /(x) is an arbitrary polynomial of degree  — 1), we obtain maximum 
degree accuracy at some fixed point x=a. There is no loss of generality in 
changing the fixed point to x=0 by the transformation x’=%x—a (and then 
dropping the prime notation). Then the desired formula is of the form: 


(2) fx P-reel®)| = YAP), 


t=1 
where P,_,.,,(x) denotes any arbitrary (” — 1+ m)* degree polynomial, and the 
x; and A; depend only on . We know that if instead of a “| , in (2) we 
b = 


were to have f ...dx, that m could be as high as . Now the problem is to find 
a 


the largest integer m for (2), for each 7. 


In the author’s thesis, [2] the impossibility of (2) was stated for just the 
first derivative and m=n, but the proof is immediately applicable to m as low 
as 2 [ie. for y=1, with m points even (n+ 1)™ degree accuracy is not possible 
for a Gaussian-type formula like (2)]. The proof is repeated here because it 
can be extended to higher order derivatives: If (2) is true for y=1, m=2, choose 


P43 (%) =(ax+5) IT (x—x,). Suppose some x;,=0. Then in P,,,(x) choose 
a=0, b=1. The right member of (2) vanishes, which implies that the coefficient 


of x in Il (x —x;) vanishes, or Py My Xq... Xj_4 Xj44---X,=0. But this is im- 
i=l j=l 

possible because there will be just a single non-vanishing term in the summation 
corresponding to7=1. If no x;=0, in P,, (x) choose a=1, b=0, and the vanish- 
ing of the right member of (2) implies that the coefficient of x in P,,,(x), which 
is simply x, x2... %,, vanishes, which is impossible. 

However, arbitrary n degree accuracy is possible with (2). Let us still 
treat for the moment the case r=1. We cannot have any x;=0, as seen by 


choosing P,(x) = [ ] (x — x;) in (2). Consider the remainder in (1) which is given 
i=1 


for any 7 in terms of the divided difference function [x x, x,... x,] by 





(3') R,(x) = te {TT «— ",) [x ty-.. mh, 
or finally (see [3]) as 

R, (2) = dele sale +i fv {[T («— *)} x 
3) x AFDC) +A)! Er (r — 1) ihe a} x 


x fF (C5)/(m + 2)! +--+ 7! IT (x — x) f° Crim +7)!, 


where each ¢; is within the least interval containing the points x and x,;. For 
v=1, (3) has two terms, and the vanishing of the coefficient of #"(¢,) at x=0, 
which is necessary and sufficient for m™ degree accuracy, is expressible by 
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n 
> 4 Xq... Xj_1 %j44--- %,=0, which for every x;=0, is equivalent to >) (4/x;) =0, 
fond 


a clit that may always be met. Thus for the first derivative, n‘* degree 
accuracy is the highest possible in (2). 

Even though for r=1, we cannot set the coefficients of both /™(¢,) and 
f"t0(C,) simultaneously equal to 0, for x=0, we may choose the points x,, 


satisfying x (1/x;)=0, which will reduce the coefficient of {+ (¢,) by reducing 


Co +, i absolute value, always subject to some restraining condition upon 
the x; to keep them within or near a given range, such as | x;| S1 after suitable 
change of variable. But we cannot expect %, x,... x, to actually attain its 
minimum yale which is here 0. For example, in a three-point formula we 


could satisfy x (1/x;)=0, |x;| 1, by setting x,=1/N, x,=—1/(N+1), %3=4, 


and also [4 Se ial =1/N(N-+1)<e for N arbitrarily large. But in practice the 
rub enters in our inability to make ¢ indefinitely small because that would 
require at least two of the points x; to become indefinitely close to each other 
and to 0. That, of course, is not possible with measurable, observable or em- 
pirically determined functions. For computable functions, we can by taking N 
arbitrarily large, x,, x, indefinitely close, and by always carrying more needed 
decimal places, reduce the computing error indefinitely, and for three or more 
points get progressively more efficient formulas for evaluating the derivative 
other than by looking for the limit of the two-point ratio formula {f (x2) — /(x,)}/ 


(%_— %). 


For r=2, (2) is shown impossible for m>2 by obtaining a contradiction by 


choosing P,,..(x)=2? J] (x—-x;) for every x;+-0 and P,,,(x)=x J] («—~,) 
i=1 i=1 
when some x;=0. For m=2, (2) applied to this same P,,,(x) shows that every 


x;+0, and, as before >) (1/x;)=0. But if (2) holds for m=2, 4 will igs 

‘ t=1 

hold for the special case of the n degree polynomial P, (x) =IT (x—x,), when 
n 


it will imply that 2, My +s Mj Kp ess Mpa Mpg +++ %y=O OF more pee 
" i+ 1,8= 
a (1/x;x;)=0. But for real x; we cannot have simultaneously p (1/x;) =0, 
t+]; M j=1 


x (1/x;x,;) =0, since the square of the former — twice the latter => (1/x?) +0. 
t+7;4,j=1 


This eliminates m=2. But m=1, or n™ degree accuracy, is anes possible 
n 


from (3), where it suffices to find x,;--0 satisfying >) (4/x;x;)= 
i+j;i,j=1 
Although again for r= 2 no higher than n‘ degree accuracy is exactly possible 
for (2), the remarks of the next-to-last paragraph for y=1 are also applicable 





* Here, for r= 2, an #; } OY also be zero, in which case the sufficient condition 
for nth degree accuracy » My... Xj %j41--- h-1%h41--» ¥n = O may be expressed 


more simply as )) (1/%;) =0. 
j+#t;7=1 
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here to indicate nae to construct formulas by judicious choice of points x;, 


always satisfying y (1/%; x;)=0 * such that in (3) the coefficients of ("+ (¢,) 
i+j;i,j= n 

and /"'*)(¢3), at x=0, will be considerably reduced, especially the )' x, x.... 
i=1 


X;-1 Xj41---*, multiplier of the "+ (¢,) which dominates the x, x5... x, multi- 
plier of the /"**(G,). 
For r>2, these same methods show at once that 1. we can always find n 
points x;, to get 2"" degree accuracy in (2), namely merely choose x, so that the 
n 
coefficient of x” in J] (x—-x,) vanishes, 2. we cannot get as high as (+r) 
i=1 


degree accuracy because (2) applied to P,, ,(x) = x’ [] (x —x;) would make some 
$=] 


x;= 0 and lead to a contradiction when (2) is applied to P, ,,_ , (x) =x7~? IT x — %x;) 
and finally, 3. even (n-+-r—1) degree accuracy in (2) ” not siontible because 


using P,,,_,(x) in (2) implies that, no x,=0, and also y (4/x;)=0, while from 


n i=1 
Fea pelt) = 2 2 TT ( (x—x;) we get also >) (4/x; x;)=0, which is impossible. 
i=l t+ 7;1,j=1 


This last discussion, on the surface, appears to leave a gap concerning the 
possibility of m-point 7‘* derivative formulas giving possibly anywhere from 
(n+ 1), (n+ 2)", ... up to (n-+-7— 2)" degree accuracy. But for any r<m those 
cases may be quickly disposed of as follows**: It suffices to show that even 


(n+ 1)" degree prrereey is impossible. This in turn is seen when (2) is applied 


to P..4(x )=*IT ( (x — x,) and P, (x)= I] («—%), which tells us that in J] ( (x — x;) 


i=1 i=1 
two cdasbiuibive ‘coefficients, not at either end of the polynomial, vanish, so that it 


is of the form (a) x*-+¢,_ 1x" 2+ +++ +¢,4,x747+0+0+¢,_.%7*+¢,_,x73+--- 

Now if («) has all real distinct roots, so must any order derivative by the well- 
known Rot.e’s theorem type of argument. But the (r—1)™ derivative has a 
‘factor of x*, implying a double root at x=0. At this point we are really through 
when we realize that any (n+ 1) degree polynomial is a special case of a higher 
degree polynomial. But this argument also applies directly up to degree (n+ 7 — 2), 
employing x’~? ]] (x —x;) and x’—-* J] (x —x,), which would imply no coefficient 

i=1 


i=1 


n 
of x? and x* in J] (x — x,), etc. 
i=1 
The final conclusion is that in Lagrangian formulas for every order of differ- 
entiation at x=0, employing the function at real distinct points x;, we can 
improve the (m — 1)" degree accuracy in Lagrangian formulas at arbitrary points 
by at most one, i.e. attain m™ degree accuracy in the theoretically exact sense. 





* See footnote, p. 217. 


** Wecan ee aneiated dispose of y = nandy = n + 1 byemploying P, (7) = I (~— *;) 


and P,.,(*)=* I (x — x,;) respectively, to obtain a contradiction in (2) for (n+ + 1)th 


degree accuracy ‘fant thus accuracy higher than the (m + 1)*® degree). 
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But even though we cannot exceed mn" degree accuracy in the exact sense, we 
can apply to (3) in general, procedures discussed above in connection with r= 1 
and y= 2, for reducing as much as is possible or convenient the terms after the 
term with /(¢,), which is dominant, emphasis being given to the term or terms 
immediately after this dominant one, since in the right member of (3) any term 
tends to dominate those to its right. For the 7" derivative, n>r-+1, it might 
be possible to take the x, fairly close together, if necessary in order to knock 
down the truncating error in finding that derivative of a computable function, 
without going nearly as far in accuracy as would be necessary in taking some r” 
difference +- h’ employing just n=r-+1 values. 

In practice, to make just the dominant term zero we might need to choose 


only one irregularly located point and thus cut down the amount of interpolation 

when /(x) is a tabulated function. For example, for y=1, to make }' (1/x;)=0, 
t=1 

X, X%q, +++) %,—, May be eament (not completely) wnaeery, and proper choice of x,, 


will suffice. When, say >) (4/x;) =0, or, equivalently, >) Kg... Xpaoy Mn++. %,=0, 


for x; fairly small the éoetiieiont Il Mj Bq sss Bi of yo (€,) will tend to be of 


t=1 
smaller order of magnitude. But we may exercise freedom in choosing more 


n 
irregularly located points x; and make ]] x, x,... x, of still lower order of magni- 
t=1 
tude, and similar ideas apply to y>1 in reducing especially the coefficient of 
fm* (Ca). 
It is worth investigating whether other choices of x; not close to each other 
(to reduce the computational error, especially for observable or experimentally 
determined functions known to only limited accuracy) will give low truncating 
errors by reducing the coefficients of {+ (¢,) and higher order derivatives. 


B. Optimal Points for Differentiation Over a Range of Arguments 


In calculating a derivative by (1) for x within any range (a, 6), (for example, 
a common range being (4%, %,), *,=%,+("—1)h, for functions tabulated at 
equal intervals of h) we might wish to minimize |R,,(x)| according to (3) in some 
suitable sense over that range, such as the least squares or Tschebyscheff (mini- 
max) sense. One difficulty is that there is more than one term, each involving 
an unknown higher order derivative. 


- Now suppose for the moment that the /“+™(¢,,.,) were known constants 
C0. Then we could consider |R,,(x)| as |co| -|S,,(x)|, the S,(x) being a linear 
combination of the fixed ratios 7,,—c,,/cy with polynomial coefficients, and a 
reasonable criterion would be to minimize S,,(x) by setting it equal to the proper 
multiple of a Tschebyscheff polynomial suitably normalized to the range. Then 
the points x; in (3) would be determinable as the roots of a polynomial x”+ a, x"—! 
+a,x"-24 ..-+ a, ,x+a, whose coefficients a, a),...,a, would be found by 
solving a linear system obtained by equating successive coefficients in S,,(x) 
with the proper multiples of the coefficients in a Tschebyscheff polynomial. Then, 
of course, the fixed points x;, although definite, would depend very much upon 
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the function by depending upon each of the ratios 7,,. But actually /“+™ (¢,,.,) 
is either unknown, or very roughly estimated, or known only to have certain 
upper bounds. The uncertainty in each of those respects causes 7,, to be equally 
uncertain, and thus the a; and finally x; cannot be determined for Tschebyscheff 
approximation which requires precise location of points. In giving up this 
approach, one is further justified by also noting that even precise values of 
{**+™ (C,,41) as a function of x could not possibly lead to a Tschebyscheff-type 
approximation to R,,(x) according to (3), because /"*™)(¢,,,,), being a function 
of x, would not determine fixed points x;. In other words, as x would vary 
over a range the change in every /“+™(¢,,.,) would alter the x; via 7,, and a,, 
and we would not be able to have a single polynomial approximation. 


Another approach to the minimization of the entire right member of (3) is 
in the least squares sense of finding the x,’s from the a,’s that yield the least 
value of the integral of the sum of the squares of the polynomial coefficients of 
the various /“+™(¢,,,,). Since the first coefficient dominates the others, this 
will show in the result. It seems feasible here also to weight the coefficients 
under the integral, possibly by some approximate upper bounds to /"*™)(¢,,.,), 
obtainable from higher order differences or divided differences (depending upon 
whether the /(x) is tabulated at uniform intervals or not). Use of fixed upper 
bounds for the entire range will, of course, give a fixed set of x,’s, and we might 
expect the weighting to be most effective and relevant in the worst cases when 
those upper bounds are approached or nearly attained. This weighted least 
squares approach then yields values of x; that also depend upon the particular 
function /(x), so that more extended computation is required in every particular 
case without the advantage of an independent location and tabulation of the 
x;'s, which we have by applying the least squares criteria to the unweighted 
polynomial coefficients of /"*™ (€,,..4). 


Before leaving either of these approaches based upon minimization of the 
entire right side of (3), we should mention another problem whose solution was 
tacitly assumed in the above discussion, but which requires special investigation 
in each particular case: In equating the right member to a Tschebyscheff poly- 
nomial or a least squares polynomial we solve real linear systems for the coefficients 


n 
a; of [J] (x—-x,). However that does not guarantee anything about either the 
i=1 
reality of all the points x;, and if real, their location within or close to the range 
of x. Meeting this last requirement might require so much extra labor and 
modification on top of the above-mentioned difficulties and uncertainties, as well 
as special calculations for each different /(x), that it seems more than justifiable 
to abandon the search for a workable criterion to minimize all of the right 
member of (3). 


There is still a workable criterion for optimization by considering the dominant 
n 


term in (3), namely - a IT («- x} {™(¢,)/n!. The extent of this dominance 
i=1 


is apparent from a fairly extensive tabulation that the writer made several 
years ago (unpublished) for y=1 and r=2, for the case of equally spaced argu- 
ments, at intervals of h, when (3) specialized into the following: 
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For r=1 
[,’ n—- n n n 
(3-1) Ry (x) =P) we p09 (G,) + LOL an fort (64), 
and for r= 2, 


_ L” (p) n—2 4(n) 2L’ (p) n—1 4(n+1) 2L(p) n t(n 
(3-2) Ry(x) = NPP wes p(t) + GE UD wt porta) (C4) + ZED pe pote (gy), 
In (3.4), (3.2) the notation for the tabular points x; followed that employed in 
Lagrangian interpolation tables, x;=%9+7h, where for the m-point formula i 
ranged bie TP 1)/2] to [m/2] (the [m] denoting the largest integer in m), 


L(p)= I[] (p14), and is the fraction of the tabular interval according to 
t= —[(n—1)/2] 


p=(x — %9)/h. Three-figure approximate values of the coefficients A a] - L(p) —_ 
L"(p) 2L'() , 14 2L(6) | re 
al? (wi)! an (n+2)1 were tabulated for p ranging from —[(m —1)/2] 
to [m/2], and inspection of the tables showed at a glance how much the first 
coefficient dominated the others. Furthermore the higher and higher powers 
of h in the other coefficients, even with higher derivatives, increase the extent 
of the dominance, since for r=1, say, h"~*/™(¢,) is followed by h" {+ (¢,) = 
he . hf (¢), and hf"*(C,) would be ~A/™(¢,), and similarly for r>1 
where the much smaller coefficients are multiplied by (after the same h"~’ for 
every, term) factors ~ 4/™), ~ A?/), etc. Thus there is ample justification for 
concentrating on the minimization of the dominant term, and by concentrating 








ary oa. . 
upon For I] (x—.,) alone we obtain criteria independent of f(x). Again two 
i=1 


approaches are possible, least squares and minimax, of which we shall consider 
only the latter in this paper* . 

The key to these optimal formulas lies in the location of the arguments x; 
at the zeros of various r“ order integrals of the Tschebyscheff polynomials, so 
that r-fold differentiation yields a polynomial identical (except for a necessary 
factor) with a Tschebyscheff polynomial of order »—yr. But the solution is far 
from straightforward because it is necessary to add polynomials of degree r — 1 
to the r-fold integral. They are far from arbitrary in that proper choice is essential 
in order to produce a particular 7” integral polynomial that will have all its 
roots x; real and distinct and also, while not necessarily within the range of 
arguments for which numerical differentiation is required, at least not too far 
outside the range. 

Although Tschebyscheff polynomials have all real roots within the interval 
of least absolute deviation from zero, it is not at all apparent whether successive 
primitives can always be found just to have real roots, even though in the study 
below it has been possible for the first seven polynomials, up to y=2. For any 
polynomial having all real roots, a necessary and sufficient condition for some 
particular r-fold primitive to have all real roots is that we can find a single 
primitive of both that polynomial and certain of its (r —s)-fold primitives (for 
every 1<s<r—1) which also have all real roots. The sufficiency is obvious, 





* The least squares approach can be developed in a way similar to that for mini- 
max, only employing instead suitable r-fold primitives of the Legendre polynomials. 
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and the necessity follows by the familiar repeated differentiation and ROLLE’s 
theorem type of argument. Establishing reality might still not be enough in 
the present case, because that might not imply anything at all concerning satis- 
factory location within or near the range. But even apart from location, it is 
easily shown that in general a polynomial with all real roots might have no 
primitive with all real roots. Thus a primitive of p(x), say P(%), may be repre- 
sented by the area under #(x) from the origin + aconstant. For p(x) of degree 
1, 2 or 3, having all roots real, from diagram a) below we see how to add a suitable 
constant K, (amount shaded) to have all roots real in the primitive; but in b) 
we see how a quartic with all real roots has no primitive with all real roots 
(dotted lines indicating root location): 











real root not possible; only 
3 roots of primitive are real. 








Thus in the subsequent work, since the existence of primitives with all real 
roots could not be assumed, the procedure for finding them had to be based 
upon observation and trial, employing special devices to reduce the degree of 
the actua] computations. 

Before entering into detail it is best to clarify the normalization employed 
and the theoretical improvement in the Tschebyscheff sense over the usual 
numerical differentiation with regular spacing in the arguments. The Tscheby- 
scheff polynomials chosen here are defined by T7,,(x)=cos(m arc cosx) and in 
the interval (—1, 1), the polynomial 7,*(x)=T,(x)/2"—* which has leading coef- 
ficient unity, deviates least from 0 in absolute value when compared with any 
other polynomial of degree and leading coefficient unity. Thus for the »-point 
formulas as x ranges over » —1 intervals of length h, or the usual Lagrangian p 
has a range of (w—1), this new # will always range between —1 and 1. From 
now on, # will denote this new Tschebyscheff polynomial variable, and the original 
x will be transformed back and forth into # in the course of the arguments. 
We seek to express R,,(x) in terms of ~ and also in such a way that the improve- 
ment factor over the usual formulas like (3.1) or (3.2) will be apparent, since we 
must take into account the difference in range of variables before making com- 
parisons. 

For x ranging from %_,—1)) tO %tq/2}, Where x,=%)+kh over (n—1) 
intervals of h, and # ranging from —1 at X= %_m_3y/9) to +1 at X= %jqjq), the 


conversion formula is x= *_ 1(,—1)/2)+ — A+ = > hp,so that(55) = (7S. 
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In terms of this new variable p where ; corresponds to x; and /(x)=q(p), 
the remainder in just LAGRANGE’S ets Beret formula is expressible as 


L I (p—#) )p™ (n)/n!, which may also be written as I (p—p,) (75 +)" h” f(t) In, 
so " that the I (x—x,)/n! or L(p)h"/n! coefficient ‘of /™(¢) must be compared 


with M I (p — B,;) (“3 y h*jn!. In other words, in allowing for the range con- 


vinden from equally spaced to Tschebyscheff points, even though IT ( (p —P;) 
[n/2] , 

can be made considerably less than L(p)= [J] (p—1), the factor of (54) 

must also be considered. Ce oA 7 


For comparison of remainders of 7‘ derivatives, apply the same procedure 
to (3) with the new # as the variable and g+™ (n,,,,) in place of /*+™)(C,,.,). 
Then in converting back to differentiation with respect to x everywhere, except 


Rod {II (¢—2,)} = P"— (p), we obtain for (3) in terms of the 


F—§ 


where we have 
new variable p 


R,,(*) = P (p) (75 i ee ” f) (C,)/n! + 7 Po) (p) (“+ are . 
(4) x ft (C.)/(m + 1)! + 7(r — 1) P&-*)(p) (5 nor 


x f+) (G5) /(m + 2)! --- +9! P(p) (">* * n)" 10* C41) /(m +7)!. 





To estimate the improvement over the usual Lagrangian methods employing 
equally spaced arguments, for which tables are already in existence, ([4], [5]), 


it suffices to compare here P”) (*) (” a " with L(p) in the dominant terms, 
then P’-)) (6) ("> = “with LY 1)(p) in the next terms, etc. 


The actual eatiniane of the upper bound for the coefficient of h*~’f™ (¢,) 
is made as follows: Since P“)(p) is to be proportional to a Tschebyscheff poly- 
nomial of degree » —r, and there is a factor of m(m—1) ... (m—7v-+1) multiply- 
ing the #"~’ from differentiating P(p)=(p — p,) ... (6—p,), we have 
(5) P()(p) = n(n — 1)... (m —7 + 1) Tit, (6), 
so that the coefficient of h"~’f™ (¢,) is m(m —1) ... (n —r+1)T,*_,(p) (Sy) 


Some idea as to how much the worst error is reduced was had by looking 
at the earlier mentioned schedules for error in equally spaced arguments for 
n=7, r=1 and for n=9, r=2. In the former L’(p)/7! was around as large as 


1/7 initially, whereas 7 - T,*(p) (7+ ty / 7!, since 7,* S 1/32, does not exceed 81/2560 
which is around 1/5 of max L’(p B)/71. In the latter, L’’(p)/9! was around 0.6 at 
the most, which was compared with 9-8 - 7;*(#) (25*)' /9!, and since T;* (p) < 1/64, 
this last product does not exceed 16/315 which is only around 1/12 of max L’’(f)/9!. 


The determination of the optimum points #; in this Tschebyscheff sense 
was carried out for the first and second derivative as far as the 8-point and 
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9-point cases respectively. The integral of a Tschebyscheff polynomial T,,(x) is 
given by the following formula*, 


(6) [Tlsasa 4 [Ba@) — Ba) sx, 


which when “anor yields 


; he ad i 2 a ee" a ‘ 

Af [Melat= Cl asieray — Gay ery + woea| tet ke 
However (6) and (7) were not particularly helpful in this task of determining 
suitable constants K,, or K, and Kg, to give a first or second primitive of T,,(p) 
having distinct, real roots ~; which were either inside or close to the interval 
(—4, 4). 

Below in Table 1 are given the first seven Tschebyscheff polynomials T,, (x), 
and the corresponding special first and second integrals found with some amount 
of trial and experimentation **. 

To reduce the labor in the computations, even though the polynomials went 
as high as the 9" degree, whenever the degree exceeded four it was found possible 


Table 1. Suitable First and Second Integrals of the Tschebyscheff Polynomials 























a ee [Tale dx+K, j (Tats) dx84+K,x+K, 
| tlh 
Oo; 1 2 2 
ME © nd 
Te :°3 6 6 
2 ss eB 4 
= — go. i: cs Sa 
2) 2%* 1 han x 6 3 3 
3 1 - 2 
a Siu nie am oe ee ; 
3 | 448—3% —3B* + 2 : -) + 4 
4+ | 824—8x*+1 <a ; +x < _ $45 3 ~~ 0-10491904 
5 | 1645—20%5+ 5% S dvtutee gx Lf = x7 — Ph 4+ 5 48 0.21695331% 
| 3 a 3 21 6 
| 32 48 4 8 3 x 
- 4 a ve = ae ee We. AP apa 
6 | 32%*°—48 74+18%?—1 - x 5 #54643 —-x 7 x 5 wot 2* 2 + 0.028 
56 7 1 8 8 14 7 
- = ae 0 an ce Oe oni a ae | 
7 | 6447—11245+ 56%3—7% | 8% 3 x°4-14% ZV+S 9* 34+ 5 x 6 +0.1444 


to suitably choose K, for y=1, or K, and K, for r=2, so that apart from a root 
x=0, the variable X=x* gave equations in X of degree no higher than four 
which had all positive roots X;, and thus the roots x; were symmetrically Jocated 
as +)/X;. Throughout the investigation of the choice of a suitable K,, or K, 
and K,, it was necessary to bear in mind getting distinct roots that were some 


* We use x and the new # interchangeably, because * is so widely used in Tscheby- 


scheff polynomial notation. 
** In this numerical work the author wishes to acknowledge the assistance of Mrs. 


GENEVIEVE M. Kimsro, who employed only a desk calculator. 
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reasonable distance apart from each other, as well as real, and if not actually 
within the interval (— 1, 1), at least not too far outside, because in almost every 
practical example we might expect to be able to calculate functional values 
slightly outside the specified range. As a rule, these primitives were not difficult 
to find, since most of the two-fold primitives did not require the employment 
of special values for K, and K,. A more painstaking investigation facilitated by 
automatic computing machinery could undoubtedly find a K,, or a K, and a 


Table 2. Optimal Points p; for n-Point Numerical Differentiation Formulas for Range 
Normalization of (—1, 1)* 














n First Derivative Second Derivative 

2 +1 +1 

3 | O, +1.224745 0, +1 

4 | +0.707107, +1.000000 +1.000000, +1.414214 

5 | 0, +0.754858, +1.047309 0, +0.831254, +1.344997 

6 | +0.464365, +0.707107, +1.076738 +0.592591, +0.943458, +1.121928 

7 | 0, 0.514311, +0.891176, +1.020436 O, +0.724242, +0.958589, +1.087005 

8 | +0.248135, +0.615639, +0.944855, | +0.263457, +0.814647, +0.983758, 
+1.000000 +1.048407 

9 0, +0.456554, +0.861866, +0.984450, 

+1.039040 


* Of course, the term “‘optimal’’ is here not employed in any unique sense, or 
most favorable in every respect. It merely means that it will yield a Tschebyscheff 
polynomial factor in the dominant term of the remainder. 


K,, that might locate aJl roots #,; entirely within (— 1,1) or much closer to the 
limit —1 or +41. The irregular points ~; of the lower degree polynomials, say 
for the 2- and 3-point cases, turn out to be regularly spaced and some have no 
advantage over the regular formulas (in fact, are identical with them), but they 
are given for the sake of completing the series*. All cubics were solved with 
the aid of a recently published table, [6] which saved considerable time. In fact, 
advantage was taken of the complete tabulation in [6] to avoid interpolation 
where it was difficult because of a nearby singularity, by choosing an exact 
argument and working backwards to the K, which thus appears as an ‘‘arbitrary 
irregular constant” in some of the primitive polynomials. All quartics were 
solved first graphically and the approximate values were refined by NEwWTOoN’s 
formula. All roots are given to 6D and are believed to be correct to within a 
couple of units in the last place: 

It is interesting to note that some of these points #; for optimization in a 
Tschebyscheff sense within a prescribed range lie actually slightly outside the 
range. This is to be contrasted with best interpolation points (Tschebyscheff 
sense) and best numerical integration points (definite integrals only, in Gaussian 
sense of highest degree accuracy), where the optimal points lie within the range. 

It is shown in numerical analysis texts that for numerical differentiation at 
the end of or outside of a given range, instead of the r+-1 terms on the right 





* However, note that y = 1, m = 3, although equally spaced, is still slightly better 
than the usual Lagrangian formula which ery yg p; at 0, +1, since the upper bound 
in the absolute value of the coefficient of h?/‘(¢,) is reduced from 1/3 to 1/4. 
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side of (3) for R,(x), we have only the first term, where @, is in the interval 
determined by x and x;. However, even though this dominant term is actually 
the only term, for calculation of the derivative outside the range (— 1, 1), the 
Tschebyscheff polynomials 7,*_,(p) mount so rapidly that these optimum points 
would be useful only very closely beyond the end-point —1 or +1 in the new 
variable p. 

One could produce a small table of first and second derivatives of Lagrangian 
interpolation coefficients, to be used in (1), based upon these irregular locations 
of p;. They might be similar to the two existing tables of derivatives of Lagran- 
gian coefficients for equal spacing of arguments, [4], [5] but perhaps with more 
emphasis upon tabulation in that part of the range where the remainder in the 
usual formula is comparatively largest. One might also tabulate briefly the 
dominant part of R,,(x), apart from the /™(¢,). Such a table might be quite 
useful as far as the first four derivatives and could go up to 11-point formulas, 
and thus in a certain sense parallel a similar tabulation by BICKLEy for equally 
spaced arguments, but which is only for x at the tabular arguments and not in- 
between. The computation could be done largely automatically as soon as the 
proper primitives have been chosen and their zeros found to a sufficiently large 
number of places. Then the coefficients in the Lagrangian polynomials would 


be obtained by synthetic division and multiplication by [J] (x;—x,)~}, there 
hi; k=l 

being a different series of m-point formulas for each 7. Finally, the r“* derivative 

is found for each polynomial in the corresponding series, and tabulated auto- 

matically for the range (— 14, 1). 

Two other possible alternative techniques for reducing the dominant part of 
the remainder should be mentioned in closing: 

1. Since the tabulation of the dominant coefficients in R,(x) for r=1 and 
r=2 for equally spaced arguments showed them to mount rapidly only near 
the ends, and since some of the above-mentioned optimal points call for some 
extrapolation * anyway, one would wish to investigate using the same Lagrangian 
formulas but for an extra end-argument outside the range, which would place 
the x further inside the new augmented range. 

2. Here we assume that we have a table of /(x) at intervals where the inter- 
polation to full accuracy requires more than m points. We first glance at brief 
tables of L) (p)/n!, for r™ derivatives, to see how to change the given argument 
p to correspond to a p’ where the dominant term is small. Then certain equally 
spaced non-tabular arguments x;+(p—>’)h will be slated to be new equally 
spaced tabular points. After we find the functional values /(x;+(p—’)h) by 
ordinary interpolation and a single extrapolation*, the desired derivative can 
be found by the usual m-point formulas or tables [4], [5] based upon equally 
spaced points. Of course, each different argument x at which a derivative is 
required would call for different new fixed points x;+ (p — p’)# and would involve 
n interpolations each time. The amount of work in those m interpolations might 
be lessened by the fact that the tabular fraction p— p’ varies only with x. 





* By “extrapolation’”’ we mean computation outside the range. But with a table 
of /(*;) the ‘‘extrapolated”’ value is found actually by interpolation, and it will thus 
be as accurate as the other values of /(;) in (1). 
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The procedure in 2. resembles that in A. for fixed point optimization where 
we chose suitable values of x;. However there is the big difference that, even 
though in A. we often might retain every x; but one as originally specified tabular 


points, the tabular arguments are still not all equally spaced and the coefficients 
d? L‘™ (x) 
A,= —aaF ~|rno? *=1>+++», must always be specially computed. In 2. the 


fixed arguments x;+ (p — p’) h are all equally spaced, so that in (1), for y=1 or 2, 


work may be saved by using existing tables of iy Lf (p) for variable # [4], [5]. 
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The Zeros of the Hankel Function 
as a Function of its Order* 
By 
WILHELM MAGNUS and LEON KOTIN 


1. Introduction 


The theory of diffraction of electromagnetic waves by a sphere requires a 
knowledge of the zeros of certain transcendental functions. 

In the simplest case of a perfectly conducting sphere of radius a, the equation 
in question is H}(ka)=0, where H? is the Hankel function of the first kind of 
order v, and k is the wave number of the incident wave. The quantity to be 
determined is the order v, with ka being regarded as a parameter. Special 
attention is paid to the case in which z= ka is real. 

We evaluate » only for sufficiently large values. For other values, we in- 
vestigate the behavior of » as a function of the argument. Specifically, we derive 
certain inequalities satisfied by the real and imaginary parts of y and the argu- 


ment z; we investigate the behavior of y as z—>0; and we obtain an expression 
dy 
for ar’ 

Although many of the results can be extended, the argument z of H} will 
be restricted for the most part to the right half plane; i.e., |argz| < =. This 
represents the physically significant case. 

As usual, a bar will denote the complex conjugate of a quantity. Because 
of the relation H(z) = H?(z), any result concerning a zero of H} can be inter- 
preted in terms of H?; e.g., for real z, the zeros of H} and H? are conjugates. 

Another immediate result is the fact that the zeros of H} are symmetric about 
the origin, since H1_,(z)=e'”*H}(z), so there is no loss of generality in letting 
the real part of the zeros be non-negative, as we shall sometimes do. 

Throughout this paper, the letters x and y will be reserved for the real and 
imaginary parts of z, and «+7 will likewise denote ». 


2. The order of growth of H}(z) with respect to v 


The preceding results, and all the subsequent ones predicated on the vanish- 
ing of H}, would be true vacuously if there were no zeros. To lend substance 
to these results, therefore, it is necessary to prove the existence of zeros. For 
this purpose we shal] apply function theoretical considerations. Specifically, 
since H} is not only analytic but entire in », we shall apply the concept of the 
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No. AF 49(638)229, and by the U.S. Army Signal Research and Development Labora- 
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order of an entire function to prove that there exists an infinite number of zeros, 
and also to obtain an infinite product representation. 


From HEINE’s formula (see [1], p. 26), for O<arg z<_a, 


co 
g(r) = = tals H} (2) = [ é#cosh vtdt. 


0 
If we let |v] =, 


co 
|g()| SS envomneter ar, 
0 


We shall show that the order of g(v) is $1; i.e., that 


le|<mMe™ forall e>0, as o>, 


where M is a suitable quantity which may depend on e¢ but is independent of o. 
It suffices to show that 

I= fi e~7ooshtt ele" di<M. 

0 

But. 

e ry (1eve 

Beg ft de+ fe“ 7m dt. 
0 e 


Since each of these integrals converges, g(v) = O(e’**) for large 9; i.e., the order 
of g(v) is S1. 

Now let »*=A. Then g(v)=g(At) =/(A) is an entire function of A and the 
order of /(A) is S}. Hence (see [2], p. 250), 


H(A) =#(0) TT (1— 3,): 


where the A, are the zeros of /(A). We note that /(0) +0 since H}(z) =-0 from 
Theorem 3.3, to be proved later. 


In terms of v, this product becomes 
1__ ,—vxi/2 py v® 
Hi = e~**4? Hh IT (1-3) 


where the », are the zeros of H}(z) for y>0. 


Expanding /(A) in powers of A and equating coefficients of the linear terms, 
we find 


+ =—/(o)/f(0), 
LEA Fou 


or 
fe gitcosht 74 
0 








» A __ gv) dv —— : 
oh 6) dane ircomtt gy 
0 
Sinee /(A) is an entire function of order <1 and is not a polynomial, it, 
together with the Hankel function, assumes each value infinitely often. In 
particular, when y>0, H}(z) has an infinite number of zeros (see [2]. ch. 8). 
Numer. Math. Bd. 2 17 
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We can extend this result to the real axis by applying the same technique 
to the function A(v) defined below. When x>0, from NiIcHOLSson’s formula 


([Z], p- 34), Q 
(2.1) h(v) = = H} (2) H2(z) = f K,(2zsinht) cosh 2v¢dt, 


where Ky is the modified Hankel function of order zero. Then if we put |»| =, 
|A(v)| S f |Ky(2zsinh 2)| e?¢* dt. 
0 
Now for any v and z for which x>0 ([3], p. 184), 


co 

(2.2) K, (z) = f e~*"* coshysds. 
0 

Therefore, 


> ro) 
|Kg(2zsinh ¢)| < f e—®*sinbtcoshs 75 < p—Bssinht  2—s*ssinht 7, 
e 0 





Thus 





—2xsinht+2Q# 
vols le [¢ oil 
2 )* sinh ¢ 
0 


To show that the order of h(v) is S1, it suffices to show that 


oo 


= [ exp [— 2 sinh ¢+ 0(2¢—9°)] di<M. 


//x sinh?¢ 





0 
Now 
of/2 


co 
sf g~Sseahs a+ [ exp [—2- sinh ¢+ (2¢)'*/*) dt. 
tt )/# sinh ¢ sinh ¢ 
0 8/2 


Since each integral converges, h(v)=0(e®**) for all positive ¢, and the order 
of h(v) is S1. 

As before, we have h(v)=h(A#) = P(A). Since cosh 2v¢ is an even function 
of v, P(A) is an entire function, and is of order <4. Then there are an infinite 
number of zeros of P(A), and hence of H}H?. If z(=<) is real and positive, 
then H? (x)= H2(x) and H}(e**x) = — H1,(x), and with our previous results we 
have 

Theorem 2.1. If OSarg z<72, there are an infinite number of zeros of H}(z). 


Since P(A) is of order $3, 
P(A) = P(o) II (1-+). 











or 


Hy (2) H(z) = Ho(2) HE (2) [TT (1 — 4) 


k 
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for x>0, where the », are now the zeros of H}H?. Moreover, equating coeffi- 
cients of v?, we have 


+t =- 20) -- 
a 2h (0) 
For z real and positive, if » is a zero of Hi, we have vf =»). Hence 


co 
2f t? Ky(2z sinh t) dt 
0 





f K,(22sinh #) dt 
0 


and 

Theorem 2.2. 
f ? K,(2x sinh?) dt 
0 





Rey pr = — 


f K,(2xsinh t) dt 
0 


To show that H} has only simple zeros, we require the following well-known 
Lemma 2.1. Let F(t) be positive, continuous, and decreasing for ¢>0. Then 


Jf F(é) sintdt>0. 
0 


Proof. 
co co -(2n+1)x (2%+2) x 
fF(@)sintdat=—> J Fi)sintdt+ f F(t) sin tad]. 
0 n=0' 2nzx (2n+1)x 
By the first mean value theorem, 
(2%+1)2 (2n+2) x 
f Fi@sintdt+ f F(t)sintdt=2F(t’) —2F(t’), 
2nx (2"+1) 2 


for some ?’, ¢’” such that 2na<t’<(2n+1)m and (2n+1)a7<t"<(2n+2)z. 
Since F(t) is strictly decreasing, F(t’) — F(t’’)>0. 

We generalize this result in 

Lemma 2.2. Let ¢(0)=0, 9’(t)>0, ¢’()=0, p(t)>0, p’(t) <0, with p'(é) con- 
tinuous, for ¢>0. Then 


Seo sin g(t) dt>0. 
0 
Let x=q(t). Since q'(t)>0, t=q74(x)=t(x). Then 


por _ f ples) © 
J PU sing(y at 7t(z)] sin xdx. 


Applying Lemma 2.1 to the function F(x) = £ completes the proof. 


Now we can prove the simplicity of the zeros. For if we let H}(z)=0 where 
x>0, from a formula of WaTSON (see [7], p. 30) 


fo] 
H?(z) oe) — Sf Ko(22sinh te?" dt. 
0 


17* 
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Since H?(z) +0, using (2.2) it suffices to show that 


co co 
Imf fexp(— 2zsinhtcoshs — 2vt)dsdt +0, 
0 0 


or 
co 


f fexp(— 2xsinhtcoshs — 2a#) sin (2ysinhécosh s + 2B) dtds +0, 
0 0 
where we changed the order of integration. As mentioned earlier, we can assume 
«=0 with no loss in generality. Since for y20 we can prove that B>0 (see 
Theorems 3.3 and 4.1), the integrand, as a function of #, satisfies the requirements 
of Lemma 2.2 and is therefore positive for each s. Consequently rae = 0 
and we have , 

Theorem 2.3. The zeros of H}(z) are simple when x>0 and y=0. 

We note that the proof is valid for the general cylindrical function of the 
form a J,(z)+6Y,(z) in which the coefficients are independent of ». 


3. Zeros of H}(z) for a fixed, complex z 
In this section we shall derive several inequalities relating the zeros » and 
the complex arguments z. For the most part, z will be restricted to the upper 
half plane. If we include the positive real axis this corresponds to the actual 
physical situation, since the wave number, which is essentially z, lies in the 
first quadrant. 
Theorem 3.1. Given #=0, there exists a positive y such that Hj,(iy)=0. 
Proof. For 6 real, we have (see [1], p. 34) 
co 
f e-teoshd HY. (i) dt = 
0 
Choose B= $ . Then the integral vanishes. Since 


— Hig (it) = Hig(— it) = Hjg(tt), 


_ 2% el? sin (B b) 
sinh B zsinhb * 





it is pure imaginary, and must vanish for some ¢> 0. 
In this case we can show that |B|>y by means of the following classical 
argument. The function H}(ze)=w«/(t)+7v(t) satisfies the equation 


7, Hi (ze) + (226?! — »2) Hi (ze) =0. 
If we multiply by H} (ze), separate the real and imaginary parts, and integrate, 
we obtain for y>0 


(3.1) f (ye — a) (u? + v2) dé=0 


bo 


where /, is a root of u(t)+-7v(¢) and 


oo 


(3.2) Fs [(a® — 92) et — (a2 — 2)] (u® +02) dt = f (u’? + v'2) dt>0. 


b 
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We now let z=». Then equation (3.1) becomes 
co 
a Bf (e®* — 1) (uw? + v2) dt=0. 
bo 


If 4420, then we must have «B=0. But then «=0 since B=y>0, and (3.2) 
becomes 


BS (A — e*') (u2 + v*) dt>0. 


This, however, is impossible if 4j=0. Therefore, 4, must be negative. 

We have thus obtained the following result. 

Theorem 3.2. If H}(av)=0 for a>0 and f>0, then a<1. 

We have as a special case 

Corollary 3.1. If Hi,(iy)=0 for y>0, then |£|>y. 

We can exploit equation (3.1) somewhat further. For instance let x>0, 
and recall that y>0. If ¢ is sufficiently large, xye?‘—aB>0. Then for the 
integral in (3.1) to vanish, since the above function is monotonic in ¢, we must 
have xye?*—aB<0. Similarly, if x<0, the inequality is reversed. Finally if 
x=0, «B=0. This establishes 

Theorem 3.3. If H}(z)=0 for y>0, then sgn(af — xy) =sgnx. 

Moreover, if y=0, then H}(x)=H?(x). But these, as linearly independent 
solutions of BESSEL’s equation, cannot vanish simultaneously. Together with 
Theorem 3.3, this proves that if y>0, then H(z) +0, a result first obtained by 
MACDONALD ([3], p. 511). 

If we choose «0 in Theorem 3.3, then x=0. 

Corollary 3.2. If His(z)=0 and y>0, then z=7y. 

If we now let x=0, then aB=0. But 8-0 since by HEINeE’s formula 

co 
Hi (iy) = =. ns (amas cosh yt dt +0 


0 
for vy real. Thus «=0. 


Corollary 3.3. The zeros of H}(iy) for y>0O are pure imaginary. 

If we now consider equation (3.2) and assume x? — y?<0, the previous reason- 
ing indicates that (x* — y*) ce?” — («2—$?)>0. This is the content of 

Theorem 3.4. If H}(z)=0 with x*<y* and y>0, then a? —f?< x? — y?. 

Further results restricting the location of the zeros can be obtained from 
HANKEL’s asymptotic formula. This expression is valid not only for the magni- 
tude 7 of z large compared to that of y and 1, when it attains greatest accuracy, 
but also for small values of z. If we consider only the first term of the expansion, 
HANKEL’s formula is 


(3.3) H}(z) = (2) Jl Fa) a4 Ry. 


We shall find conditions under which the remainder R has magnitude less than 
unity, whence H}(z)=+0. First we shall use WEBER’s form for the remainder. 
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The derivation of this formula for real values of » is reproduced in WATSON 
((3], pp. 211—212). This can be extended to complex values by only slight 
changes. Without going through this lengthy analysis, we merely state the 
following inequality for R. If 2|z]|=2r2a—4>0, B=0, and x>0, then 














2 1//,_ #—4)-#4-1 | F(a+4) [2 
*—a\-Sy lpeah 
But it can easily be shown that when («+ 4) (a —4)<r, 
a 2r—a+} sinh B x 
[R| 23 > * 2r—(a+1)(2a—1) Bau ~* 
This implies 
Theorem 3.5. If H}(z)=0 for «>}, B=0, and x>0, then 
2|z| <(«+4) (2 4 wee. 








We can complement this result insofar as « is concerned by considering 
another form for the remainder due to ScHLAFLI [(3], p. 219). When |a|<# 


and — — = <argz< 3m, 


COS ¥ 7 
COS @ 7 


Ja?—4\c < cosh Bx Ja?— 2I¢ 
2r = |cos a 2| 2r 


(3.4) |R| Ss 


where c=1 when y2O and c=|secargz| when y<O. Then if H}(z)=0, we 
must have |R|=1, which gives the following result. 


Theorem 3.6. If H}(z)=0 for |a|< 3 and — — ~ <argz< 4 ——, then 


HE > 
4 











or< c cosh B x 











where c is defined as above. 
Since |«|<#%, we obtain the following less accurate but simpler inequality 


for y=0: 
cosh B x 
= |cosa 2| © 


4. Roots of H(z) = 0 for real, positive values of z 


Most of the results of the preceding section are valid only when y>0. In 
this section we shall restrict z to the positive real axis. 


Lemma 4.1. If H}(x)=0 and x>0, then 


1 4 = 
fis oes = 


Proof. We consider the self adjoint form of BEssEL’s differential equation 
for H}(x) and H2(x). After multiplying each equation by the other function, 
subtracting, and integrating, we obtain 








[AF Ro moa=i[mo 7 HO - HO 7 Hol? 


x 
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or 


4iaB f \Hr (ol? am 4 ps 


from the behavior of the Hankel functions at infinity, and the fact that 
H} (t) =H; (2). 

Since the last integral converges, we have 

Theorem 4.1. If H}(x)=0 for x>0, then «B>0. 


As far as « is concerned, this tells us only that | «| >0. We can improve this 
inequality, but first we require the following well-known extension of Lemma 2.1. 


co 
Lemma 4.2. Let f /(é) cos bt dt exist for b +0, where /(#) is twice differenti- 
0 fo) 
able for ¢>0 with /'(#)<0 and /’(t)>0. Then f f(t) cos bt dt>0. 
0 
Proof. Integrating by parts, we obtain 


fie cos btdt = + f(t) sin bt|>— s fre sin btdt. 
0 


0 
Since f /(t)cosbtdt exists with /’(t)<0, then both lim f ()=0 and f(t) =O(¢) 
0 oo 


as t->0. Then é 
fi cobiat=— 5 fr ) sin btdt, 
0 


and an application of Lemma 2.1 aieniien the proof. 
We are now able to prove 
Theorem 4.2. Let H}(x)=0 with x«>0. Then |a|>x— )/x/2. 
Proof. From NICHOLSON’s formula (2.1) for x>0, 


Re e H}(x) H5(x)| ms fre cos 2B tdt 


where 
f(t) =K,(2% sinh t) cosh 2«¢ 


= f exp[- 2xsinhtcosh s] cosh 2atds. 
0 
After differentiating twice it can be shown that the conditions of the preceding 
lemma are satisfied when a>0 and 4x2>4a%+1+)/8a+1, whence H® (x) 0. 
The obverse then implies the theorem. 


If we restrict |«| to less than $, we obtain another inequality. 

Theorem 4.3. Let H}(x)=0 with x>0 and |a|<}. Then x<ftanaa+}. 

Proof. Using a formula due to Watson ([3], p. 445), we find for x>0 and 
|o.| <<} that 


stm zen f Kral (2x sinh 2) e e~2iBt qt — H) - 5 (x oe x) — H} (x )H ~+(4) 
= [HS (fhe — [He [Bem 
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Then 


co 
(4.1) cosa [ Kyq(2%sinh t) cos 2B tdt= = (| A? |2e°* + | Hp |? oF) , 
0 


co 
(4.2) sin am { Kyq(2xsinh t) sin2Btdt= = (||? e6* — |H} |? e—F) . 
0 
After subtracting we obtain 
P co 
[Hi (x)[*e-P* = f K,_(2%sinh #) cos(2B t+ a2) dt. 
0 


After dividing both sides by cos «x, integrating by parts, and applying equation 
(2.2) for K,, we find that when x>ftanaz+3, Lemma 2.1 is applicable. 


We can close the « interval under consideration in the following manner. 
If we let yx=4+7£, then 


a = 1(H} H3_, — H}_, H?) =i(H} H*, — H1, H?) =|H}|? e~?* + |H3|2 6, 


where the argument of both Hankel functions is ¢. This is as easily proved for 
a= —4, and incidentally provides an extension of (4.1). Now if we divide both 
sides by ¢ and integrate from x to oo, applying Lemma 4.1 we find 


00 
aaa te flor, 
whence r 
Theorem 4.4. Let H!(x)=0 with x>0, a=+ = Then |p|> =. 
If we apply Lemma 2.1 to equation (4.2), we obtain 
Theorem 4.5. If x>0 and |a|<4}, then 


sgn [|H3(x)| e°* — |H¥(x)[] =sgnaB. 


5. An application of HURWITZ’s Theorem 


We can extend the results of the preceding section by using the following 
theorem due to Hurwitz ([2], p. 119). 


Theorem 5.1. Let /,(z) be a sequence of analytic functions which converge 
uniformly to /(z) =0. Then 2 is & zero of /(z) if and only if it is a limit point 
of the set of zeros of /,,(z). 

In particular, if /,(z) is the partial sum of a power series, we can use the 
following result due to CaucHy to get a bound on the zeros of the limit. 


Lemma 5.1. Let /,,(z) =a@)+4,2+ 4,27+---+a,2"=0 with a,+0. Then 


1 -1 
|z|=> [1+ Te] ™* |4,I| ; 
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Now for x>0, 
- Hy (x) H3(2) )= [Ke (2x sinh t) cosh 2¢tdt 
ane be a, y?* 
n=0 
where 
co 


a, = Gayl ij (22)2" K,(2xsinh?#) dt 


after reversing the order of octiniies and integration. Applying Hurwitz’s 
Theorem, we obtain 


Lemma 5.2. If H}(x)=0 and x>0, then 
|y| => 1+ 2 maxa, mJ 
0 n" 
We note that 


a, < f Ko(2x sinh #) cosh 2¢ dt 
0 


2 
= = |Ai(*)|? 
and 
2 
a= = lo (2)/?. 


This results in the following simpler inequality. 
Theorem 5.2. If H}(x)=0 and x>0, then 


Hi 
bl >[1+| ert | 


3 \24-4 
2 (1+ =) 
I> > i ee ma 








which implies 





and, if x>4, 
8 2)-4 
io] > [1+ [etary 


The latter inequalities are obtained from an application of the asymptotic 
formula (3.3) and (3.4). 
These inequalities can be improved at the expense of the range of x. 
Theorem 5.3. Let H}(x)=o. If “ee then 


|»| 2 (1+ we 





u% aaHIer 
Proof. It can be shown that arc siah s<s* where k=; and s=0. Then 


oo 


<2 [mn (2xs)ds 
(2)! ° res. 
0 








: 
_ 
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According to [3], p. 388, this is 
22-2 


= —?__J*(kn++)= 
(2m)! x2*"t1 ( ii 2) 

Forming 6,.,,/b, we find from the fact that I"(x) increases for x21.5 that 
bn41<0, for n=1 if x**>0.5. For n=0, b,>by if x7*>0.606 ... =2I?(k+4)/z. 
Thus a,, <b, for all m and the theorem follows from Lemma 5.2. 


The resulting inequality is stronger, as asserted; for when x>0, 


co 
Hi (x)? = 5 f K,(2xsinh?) cosh 2¢ dt 
0 


> = [ Ko(2xsinhs) cosh t dt 


6. Behavior of the zeros as x > 0 


The inequalities | «| > x —/x/2 and «<< tan a 2+} which we found previously 
give us no information when 0<*<}. We can explore one end of this range 
by letting x tend to zero, which enables us to employ approximation formulas 
of a relatively simple nature. 

First we — that «—>0 with x. For if we multiply both sides of Rove: S 
equation by Hi(#), integrate the real part, and apply Lemma 4.1, we obtain 


(6.1) flap ae atm Ae 


x 


We anticipate the results of the sequel by noting that # is an increasing function 
of x (Corollary 7.1) and thus remains bounded as x->0. So does «, since « and 
B can approach infinity only simultaneously (Theorem 8.1). Now consider the 
integrand for small ¢. An examination of the series expansion for H}(#) and its 
derivative reveals that the latter dominates the integrand, causing the integral 
to diverge to — co as x->0. Since «B>0, we must conclude from (6.1) that «0. 


With this we can prove the following. 
Theorem 6.1. Let H}(x)=0. Then 


lim v1 « log x = 0. 
270+ 





Proof. 
(6.2) H, (x ‘- ls sin vm i a) 
= AZ) Te +e (SZ) T(—»} (14-.0(2) 
as long as v is bounded away from —1, —2, —3,.... Hence 


lim {e-a¥6 T"() — exp a log = +Ba| \r'(- »|} =0. 
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Sincea—>0, |I"(v)| ~|I"(—»)| and we find that lim (alog = oa Ea) =0. It remains 


to show that 8-0, for then r')~- — and the proof would be complete. This is 
shown in 

Theorem 6.2. Let H}(x)=0 with x>0+. Then »->0. Moreover for some 
integer k, B logx+ka=o(vx). 

Proof. We know that «—>0. Suppose B+>0. Then since f decreases with x 
(Corollary 7.1), 8 approaches a limit 8, which we may assume positive. Applying 


(5) ~ el? from the proof of the preceding theorem, to (6.2), we obtain 
ePol2 T"(i Bo) [cos (Bo log 5 — isin (6o log 5) + 
+ efor? T"(— 4 By) [cos (Bo log 4 + isin (Bo log 5) =0(zx). 


Choose a sequence of x’s converging to zero such that sin (60 log 3) =0 for each x. 
Then J"(¢B,)+I°(—7B,)=0. Similarly choose another sequence of x’s such that 
cos (60 log 5} =0. Then J"(18)) —I"(— 1B )=0. These results are incompatible; 
thus By=0. 

Therefore from Theorem 6.1 we have (=)"" —1. Also I"(v)~ x Then 
from (6.2) 


ow a Oh sa x 
o(x) = rs [cos (8 log 4 7 sin (8 log 5) 
— oe “i (8 log 4 + isin (8 log 5) 
= — sin (6 log — 3) 
Say ame! Som. fase BP - Zz: 
: * (+) (tae + Blog 4 
for some integer k, depending on which zero » we take. Then # logx+ka=o(vh). 


7. The derivative of 
In this section, we shall derive a formula which was first obtained by WATSON 
(see [3], p. 508), 
oa ref K (2c sinh #) e~*** dt, 


where v is real and c is a positive zero (in z) of the cylindrical function 
J,(z) cos s+ Y,(z) sins for constant s. This form immediately excludes H}. We 
shall generalize it to apply to complex values of z and », as well as to any 
function of the form C,(z)=c, J,(z)+c,Y,(z), where c, and c, are independent of 
both z and ». 

All the results of this section will be corollaries of 

Theorem 7.1. If C,(z) is of the above form, with x>0, then any zero » satisfies 


(7.1) 2G [Ke (2zsinh t) e~*"*dt=1. 
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Proof. Let B,(z)=, J,(z)+5,Y,(z) be a function of the same type as C,(z) 
and linearly independent of it. Then 


6 = b,c, — bac, +0. 
The Wronskians of these functions with respect to v and z are respectively 


ac, OB, _ ey, aj 
B,> C, ie (15 —Y, =) 








co 
a‘ 4° [Ko (2zsinh t) e~2*# dé 
0 


4 


([2], p. 30), and 
p 4 _¢ 4B, _ 28 


» dz " dz nz 
If C,(z) =0, then 
dC, dv aC, _ 
Bote + ae oe = 
whence S +0. Each term can now be replaced by one of the above Wronskians, 


giving us the theorem. 

We shall restrict C, to H}, although many of the following results are valid 
for the general cylindrical function. 

Corollary 7.1. Let H}(x)=0 for « and x positive. Then «’(x)>0 when a=} 
or x=; and f’(x)>0. 

Proof. Solving (7.1) for & and separating real and imaginary parts, we find 








(7.2) bod =2% ia * [ Ko(2sinh # e~**4 cos 2B t dt 
0 

and 

(7.3) oF =2x ci * [ Ko(2xsinh 2) e~ 24! sin 2B tat. 
0 








If we express the kernel as 
f(t) = Ky(2 sinh t) e~ 2*# 


oo 
v4 = fexp[— 2xsinhtcoshs — 2aé]ds, 
0 


the corollary thus follows from an application of Lemmas 4.2 and 2.1 respectively. 
Moreover, with «>0,*we have from (7.1) and (2.1) 


1S 2x|r'| [ Ko(2xsinht) dt = x|>'| |H3(x)/* 
0 


Since 


2 \s 1 
\He(2)1 $ (5) (1+ gz) 
from (3.3) and (3.4), we obtain 
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Corollary 7.2. If H}(x) = 0 with x>0, then 


|»’(x)| = EGF , 
which implies 
|»’(x)| 2 as 
a ( 1+ #) 
We note from the first inequality that lim | »’(x)| = 00. 
A similar inequality is obtained from (7.2) and (7.3). For from these, 


co 
|a'| <2x|>'|®f Ky (2xt) e~2#¢ dt 
0 
= P(S — 1) tare cosh” if a>x 
x x 
= |>"l*(1 — =) tare cos—, if a<x 


(see [3], p. 388). The same inequality is satisfied by |A’], and we obtain the 
next result, where y=arc cos (h) = , the hyperbolic function being used if | «| > x. 
Corollary 7.3. If H}(x)=0 with x>0, then 


ay 
ax 


sin (h) y 


<= 2y 











If we consider the k-th zero », and let k-> oo, since a,—> oo (see section 8) 
we find that lim |%| = 00. 
—> 00 
These conditions on |»’| give us no picture of the trace of the zeros in the 
y-plane as x increases. We can remedy this by considering f as a function of « 


and determining the behavior of . 
Corollary 7.4. Let H}(x)=0 with x>0 and «>4. Then 
dp 2B 


da = V8a—1 ; 
Proof. From (7.2) to (7.4) 
ap us {tl sin 2B tdt 
da 








Tilt) cos 2B t dt 
6 


where «>. After integrating by parts, we find that sd <286 whenever 6 
satisfies 


fro + 6f'(t)]sin 2Bidt<0. 
0 


Now if 6>——— , the kernel is negative. Its derivative is positive if 6> ———_. 
2(e-Fa) . Sa—1 


We note that this value exceeds the preceding one, and applying Lemma 2.1 we 
find that < <2£6. This implies the corollary. 
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These results have so far been restricted to real values x. For complex values z, 
we get the following analogue of Corollary 7.1. 


Corollary 7.5. Let H} : =0 with x, y, «, B positive. Then 


op 1 
* 5p Gq >0 when 27: 
and 


x 2 + y2 _ ~ > 0. 
Proof. From (2.2) for «>0 


co oo CO 
f Ko(2zsinh t) e~?"*dt = f f exp[2zsinhtcoshs — 2vt]dsdt = ffe- —*4dsdt 
0 00 


00 
where 
p(t) =2xsinhtcoshs + 2at 
and 
q(t) = 2ysinhtcoshs + 2Bt. 


Substituting this in (7.1) and separating real and imaginary parts, we get 


co co 
XO, — YB, =2|z'|*f fe-? cosqdtds 
00 
and 
co co 
xB, + ya, =2|z0'|?f fe-Psingdtds 
00 


where we changed the order of integration. An application of Lemmas 4.2 and 
2.4 to the integrands then yields the result. 

These inequalities can be expressed more simply in terms of the partial 
derivatives with respect to the polar coordinates y and of the z-plane, if we 
make use of the Cauchy-Riemann equations for the we ee y(z). 


Corollary 7.6. .% H? (z) =0 with x, y, a, B positive. Then +e <0 and — so >0. 


If in addition « 2+ —, then oa >0 and x >0. 


8. Behavior of the large zeros 
In this section we shall approximate the large zeros by considering the be- 
havior of H}(z) for |v| large compared to 4 and |z]. 
The first Hankel function may be defined by 
H1 (2) = Jao) le) 


tsinv x 





It follows that 
— ivsiny 2 H}(2z) = A(v) + D(») 
where 
- ion 


At) = riz “5 + To) 





and 
eM = D,(v) + D,(v). 


——— 
- Tovey + Fete 





Di) =»>'4 


k=l 
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Using STIRLING’s approximation for the gamma function, we can approximate 
the large zeros of A(v). Then we can show that |A(y)|>|D(»)| on a circle of 


radius fiog J ~* centered at each zero. By RovucuE’s Theorem, the zeros of 
A(v) approximate those of H}(2z). These results are summarized in 


Theorem 8.1. Let z=re’? and 0<argy<z. If |»| is sufficiently large, the 
real and imaginary parts of the zeros of H(z) are given by 


cn a(S e+ Son 2) tere 


B=n(n +3) [log {#482 )-* [1+ 6,] 





where » is a large integer and e,, 6,0 | log a7) , 


log n 


This shows that the distance of two consecutive zeros of H}(x) tends to zero, 
and that the argument of the zeros tends towards 2/2, although their real part 
tends towards infinity. Therefore the behavior of the zeros as derived from 
the well known formulas of VAN DER Pot and BREMMER is not the final one. 
According to these, we have for large x (i.e. for x >1), the approximate expression 


¥,~x+a4r, 
for the n-th zero », of H}(x), where for large n 
T, = $[32(n + @)]b em" 
and where 
To = 1 .856 ef 7/8 
T, = 3.245 ef 78 
and so on. 
9. Survey of results on the zeros of H?(x), x>0 


Let x be a fixed positive number. Let y=a-+78 be a complex number and 
let H} (x) be the Hankel function of the first kind and of order v with argument x. 
We have: 

(1) The function H} (x) is an entire function of »y which vanishes-for infinitely 
many values of y= +¥»,, k=1, 2, 3,... 

fo] 
(2) era HY (x) = HB(x) [] (1-22), 
A 
where the », depend on x. 

(3) For values |v| which are sufficiently large compared to 1 and to x, the 

zeros of H}(x) are given by », = +(«,+78,) where 


=F (r+ a)feele+ a) alt [+ oC ee") 
peal foros) Sh OCR 


and where » is a (large) positive integer. 
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(4) For the zeros », of H}(x), with k=4, 2,3,..., we have 





/ 
- ft K,(2xsinht) dt 
1 
Rea = — #2 
= { Ko(2* sinh?) dt 


where K, denotes the modified Hankel function of order zero. 

(5) The following inequalities hold: If y=a+<«f and H®(av)=0, where a 
is real and positive and B>0, then a<1. 

If y=a+72B and H}(x)=0 where x>0, then 


a) aB>0O, 
b) al >*— 2/2, 
c) xSPtanaa+4 when [a] <}, 


1 
d) |A| >> when le|=—, 








e) lo] > [1+ ae nea 
f) lim »=0. 


(6) If y=a+7B and H@)(x)=0, then a, B may be considered as functions 
of the positive real variable x. We have 
1 


da 1 
az 9 whenever a2— or %az- 


SE >0 whenever a>0O and x>0. 


dx dx P 





z—>0 





ap — 2p when a>. 


da = VW8a—1 8 
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Approximation to the solution of partial 
differential-equations by the solutions 
of ordinary differential-equations* 

By 
RUDOLF F. ALBRECHT 


I. Introduction 


This paper is concerned with general methods to approximate the solution 
of a given boundary-initial-value problem for a linear partial differential equation 
in two variables x, ¢ by the solution of a system of ordinary differential equations 
with respect to ¢. Such an approach is needed to obtain approximate solutions 
of partial differential equations by means of an analog computer for such com- 
puters can only solve ordinary differential equations. 


More precisely, the problems we will consider are of the following type. Let 
@,(z,t), k=41,2,...,m, 
b, (x, t), i= 0,1, ..<,%, 


and /(x,¢) be real and continuous functions of x and ¢ on the rectangle R: 
Os*Sa, 0StsSd, with 


A,(x,t)=1, 5,(x,t) +0 on R, 


and let 
= ah 
T= Da, (%,t) = (I, 4) 
k=1 
and 
. al 
X = 914, (x,t) on (I, 2) 
l=0 


Then we will consider the linear differential equation 
(T+ X)u=f (x,t) (I, 3) 


and assume that it has in R: 0<*<a, 0<t<b a uniquely determined solution 
u=u(x,t) which satisfies for 0< «<a the initial conditions 


2 u(x, 0-+) = 8 (x), k=0,1,...,m—1,** (I, 4) 
and for 0<¢<b the boundary conditions 
C,u(x,t)=c,(t), v=1,2,...,". (I, 5) 
* This work was supported by NASA Contract No. NAw-6557. 


we 3P u(x,O0+) = u(x, 0+). 
Numer. Math. Bd. 2 18 
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The functions “§,(x) are real and continuous in 0<%<a, the operators C, 
are linearly independent and defined by 


™ , 
C= Dieailt) sail eno for v=1,2,...,S, 
I, 6) 
ly 9 ( 
C= 2 eld laa. for y=s+4,...,%, 


with /,S"—1, and the functions c,,;(¢) and c,(¢) are real and continuous on 
ostsb. If 
O= Xp <M OA yy < ty =a 


is a given equidistant partition of the interval [0,4], we substitute for (I, 3) 
finite-difference expressions with respect to x at the points x; such that the 
solution U;(¢) of the resulting system of ordinary linear differential equations 
with respect to ¢ is an approximation to the functions (x,, #). 

In the following section we will describe general methods for finding such 
finite-difference expressions. In section III we will use them for the approxi- 
mation to the solution of the partial differential equation. For a restricted but 
important class of problems we will give an estimate of the error and a con- 
vergence test for our approximations. In the last section we will consider some 
applications to illustrate the principle of our approach. 

There are already several papers written on particular cases of the methods 
described here (see, for example, references [1], [2], [3], [4]) but, as far as the 
author knows, no general treatment was given and no attempt has been made 
to prove the convergence. 

Finally it should be pointed out that the methods described here may be 
extended to more independent variables and to systems of partial differential 
equations. A variant of our approach with a complete convergence and existence 
proof was given by E. RoTHE for the heat equation with a nonlinear source 
term [5]?. 

II. Remainder Operators 


Let {«, B, y, ...} be the set of elements «, B, y, ... represented by real func- 
tions a(t), B(t), y(é),... of ¢ which are continuous on 0S#Sb. This set is a 
commutative ring with respect to ordinary addition and multiplication of the 
functions. In particular, the field of real numbers is contained in it. Using 
this field as a domain of multipliers and introducing the norm 


|| «|| = max {| (¢)|} 


#€ [0,5] 

for any «, we obtain a Banach space ® with the elements «, B, y, .... 
Further, let {f, g, h, ...} be the set of elements /, g, 4, ... represented by real 
functions f(x, t), g(x, #), h(x, t),... of x and ¢ which are continuous on R. This 


set is a commutative group with respect to ordinary addition of the functions. 
Using ® as a domain of multipliers with respect to ordinary multiplication of 





1 Substantially the same method is described (but not proved) in reference [6]. 





Approximation to the solution of partial differential-equations 247 


the functions and introducing 
= max t 
Il = max {1/(2. 9) 
for any element /, we see that ||/|| satisfies all the norm axioms except that 


l= fll S [lel] [NAll- 
Using ||/—g|| as the distance between any two elements /, g we have a metric 
space 8 with the elements /, g, h, .... 

The convergence in both ® and % of the functions is uniform and the spaces 
® and B are complete. 

By %, we denote the subspace of % consisting of all functions which have 
continuous partial derivatives with respect to x through order 9 on R. We 
define 8,=%. Finally, let h be a real variable, 0<hSa'<a, a’ fixed. 

Then we consider a class of operators depending on a point x;€ [0, a], on h, 
and on the non-negative integer 9, which we term ‘‘remainder operator of degree 
o at x=x;,” and symbolize them by P?, Q{,,... or simply by P®, Qf,.... We 
define them by the following properties. 

P 1: they are defined on a non-empty subset of 8, 

P 2: for any / any /€ of this subset, and for a fixed x; 


PAH y(Pistie.t), viPGELOER, 
P 3: for any ¢, OSS), 
lim y ( P£;f;h,t)=0 


For any «€ 8, for the same h and any two different remainder operators 
Pf, Qf, defined on the same subset of 8, we set 


(a PR) f=a(PPf), (PP+ O)t=Fei+ Of. 
Then a P® and P®+ Qf are remainder operators also. 


An sateaaieah class of remainder operators involving a given differential 
operator Y and defined on 8,, g=, can be constructed as follows. Let 


Y= => ¢ (x, t) | ox 4 
i=0 
and 


C, => [Fy Lsn—1; v=1,2,....0; wn; €€[0,a] 
be given operators with domain %,. We assume 
g.€B, 1=0,1,...,"; g,(%,t) +0 on R; 
C,7ER; v=1,2,...,.@; 7=0,1,...,45 C1, (4) 0 on [0,8]. 


For any fixed x;€[0, a], any %;,,=%;+h, m an integer, x;,,€[0,@], and any 
/€%,, we denote 


hee t) by hitus 
ZrHaet) by Ma A=0,4,-:0 Mam hen 
81(% itp’) by = itp? 


Date +p t 





oat Sl. ten by Yieul- 
18* 








248 Rupotr F. ALBRECHT: 


Further, let «;,,, B;+,, and y, be ER, and let o,, 1,02, tT, be non-negative 
integers such that %;_4,, %j42,, Xia,» %i+2, are €[0, a]. 
Then, for any {€%,, we consider linear expressions of the form 


‘ ~ 0, or 
K" 0544 Yj italitn o II, 1 
sy, te tuft 2 Bis fis 3 Dy, C, 4, if Mira. 
x=1 


Using Taylor polynomials with respect to x at x=x; and remainder terms, 
expression (II, 1) can be written 


> h" City de, soe oe “iy — fern 4 »} Bisa > lal f+ 


u=—0O, u=—O, y=0 





oe (IT, 2) 


S iy, > 5 {ery fie) + Rif. 
=1 7=0 v=0 

We want to find quantities a,;,,, B;,, and y, such that the coefficients of 

t;, fi, ---,# in (II, 2) are identically zero. In addition we assume that our 

problem i is formulated for the minimum number of quantities «,,,,, Bji_4, 7, and 

that «;=1. In other words, we have to find a solution of the linear system of 

o@+1 equations which one writes in matrix form 








¥ whe EWE H” 89 itu ‘eee wh eave rork h* c,.4 eoneee a d : 7 TA” go i | 
: : : City : 
2 “(uh (ah) s. (oh) dT: 
-h De. itp eM eee 7 WHY Out mT “(al : h” g. 3 
b i= *: fo : (II, 3) 
Bi+n : 

















where in the (0 +1) X(@+1) matrix of (II, 3) elements of the first and the t+ 1st 
line are written, OS tSo, and where we use the convention 


&i+p=0 for t>n, 
and 
C.,=0 for t>/,. 
The rank of system (II, 3) has to be 


eti=ntntatatit{, 


From the determinant of the (9+ 41) x(o+1) matrix in (II, 3) we can factor h* 
+1) 


out of the t+ 1st line, r=0,1,..., 0. Omitting the factor h we have a 
determinant which is a polynomial in /. It is different from zero for all suffi- 
ciently small values of h, if it has a nonzero constant term. By assumption all 
the g, ;+, and c,,, are +0. After factoring them out we have the following 
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determinant 
is, Pie's die © Bésabcsch. 6: dis 
0 
1 
; F ; o 
Het HV 0 : ° (II, 4) 
eivues 6.204571 5. 200.5 oo 
. n—tb, 
. 
i me aol 
(e—n)! e! (o—4,)! 


which is now independent of the coefficients g,(x, ¢) and c,;(¢). If this determinant 
(II, 4) is different from zero, then uniquely determined quantities «;4,, Bi4.an, YnCR 
will exist. For applying CRAMER’Ss rule we see that substituting the right hand 
column of (II, 3) for any column of the determinant of the coefficients we can 
+1 

again factor out h . As a consequence these solutions of system (II, 3) are 
uniformly bounded with respect to A and are continuous functions of ¢ on [0, 3], 
and they are not all identically zero because g, ; = 0. 

A necessary condition that the determinant (II, 5) be different from zero 
is that 


%1+%Se+1i-n, 


and that all the /,’s are distinct. 
Using these solutions «;,,,8;,,, and y, of (II, 3), expression (II, 2) can be 
written 


Ri f= ne a City pm Siit+u hr! ae a a Bitntitn +0, oF 
p=—-O, 1=0 au=—O; (II, 5) 


wo Lx ’ 
+ p Vx > Cxj hx e(”) , 
x=1 j=0 


In the sum inside the brackets the coefficient of any of the e’s is uniformly 
bounded with respect to 4 and continuous on [0, ]. The remainder terms ¢;,,, ;, 
€4,, and e® are ER, and 


lim e¢; = lim e,,.=lim e =0. 
po *t4! pig 't+8 pg? 


Hence 
lim h~¢ R3f =. 


In addition, (II, 5) is of form P 2 and is defined on 8,. Thus expression (II, 1) 
with these coefficients «;+,,,8;;,,7, defines a remainder operator Rf. 

The method just described has the disadvantage that we have to solve a 
linear algebraic system for each particular operator Y. Under somewhat more 
restrictive assumptions there is another convenient procedure to find remainder 
operators of the desired type. This method is based on an idea of H. SASSENFELD 
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[7] and originally used by him to approximate the solution of ordinary dif- 
ferential equations. 
We assume that the functions g,(x, ¢) of the operator Y are €%,, and that 


a ' 
Cash 7=0,1,....@—1. 


Then let {€ %,, be arbitrary and set 
nn 2 n a 
vim ue Step a6 


The functions ¢,(x,¢), e,€B,, are certain sums of the coefficients g,(x,¢) and 
their partial derivatives with respect to x. They can be found by integrating 


x & * 
f J ai J Yfdx" by parts until no derivatives of { appear under integral signs. 
0 
n-times 
If we differentiate the result m-times with respect to x we can obtain the form 
(II, 6). 
For an arbitrary {€%,, we consider the following expressions with remainder 
operators “/) P? or “ P? respectively. 
(a) if the operators C;,, are involved: 


; 1 T; 2 0 if kj 
(*1) pe _— h* Cie (*) (i) itali + , , II, 7 
ee” ent ene epee e 
(b) otherwise: 
O) Pe f= p hk al | p a Bisalisn: (II, 8) 
b=—O, See 


In both cases k=1, 2,...,. 0, and t, depend on k and 7 in case (a) and on 
k only in case (b); 0, and t, depend on 7 in case (a). All coefficients «;,,, and 
8:4, are constants. 
Using in case (a) the operators “A;, 9A;, “IC, defined by 
OA t= > MMe hs, for R24, 


a=—0O, 


4, f=— > OB: ahi+n 


oncy =| 0 if kSj, 
—WCif if k>j, 


and in case (b) the operators A;, A;, defined by 


MA t= > Me alien for k=1, 


B= 0; 
A; f= — 2 Bits hitn 
for an arbitrary /€ 8, equations (II, 7) and (II, 8) become 
Cpe s— NA. f) NA ¢ — FNC fF, 
pes — 4.) _ 4 f. 
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The operators “A; are assumed to be the same for all k’s when j is fixed while 
the operators “A; and “?C (k=1) depend on k; also the operator A; is the 
same for any k. Of course, these are restrictions but they enable us to find 
very easily remainder operators involving the operator Y. We have already 
mentioned that remainder operators of the same degree, at the same point x,, 
with the same value of 4 and the same domain of definition, form an Rt-module. 
Thus, for /=0, 1, ..., ”, we conclude from 


(1) pe f — 094, — 4 f — MAIC f 
and 
(n—11) pe f ps (m—bI A, fn) oe 4 f we ("LAC f 


by subtraction because of 


IC="--NC for n>l+j 
that in case (a) Q9, defined by 
CAO f = iA, fe") _ (hi) 4 s9—D __ 0 if m>l+]7, 


Ct if n<l4j, 
is a remainder operator, and similarly that in case (b) "Q$/, defined by 
Qe f = A, j sie wm, fim, 


is a remainder operator. Let us assume that the operators “Pe, “Pe are of 
the type described in our preceding method. Then by (II, 5) for any polynomial 
f in x over Rt of degree So 

(ki) pe f per (*) pe f =0. 
As a consequence, "0? and “GQ? defined by 


0 if m>l+7, 


("AC fiat if n<l +4, (II, 9) 


(NOE f _ 4, fo oi (LA f a 


and 
f= A, /9— 0A, (II, 10) 


are remainder operators for any /€ B,_,41- 

Finally, we show how the remainder operators given by (II, 1) and (II, 9, 10) 
apply to the partial differential equation (I, 3) with the boundary operators (I, 6). 
For Y=X, for u(x, t)€B,, and for =0 or =a, respectively, we have 


0, or 
Ss 
t s Wny C,, OF 
Peu= » Wa; 4.4(— Tia u@ + hiu) + pa Bit-n“itnt 2, 
ei | — " 
Dd WY Cys 
x=s+1 


with 
Ti+n sme [T@] ron; 5° 


Using the operators of (II, 9, 10) we have 


* » (NOG (e, 4) = “PA, (Tu — f) + » (LIA (eu) + D “MC(eu)"r, 
i=0 l=n—j 


i=0 m7 
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or 
—Z O(c.) =A (Tu) + FPA (en). 
If e; € Bo nth then by 
SOR f) and FOE) 
1=0 l=0 


remainder operators for {€B, are given. If m=2mn’ and if X has the self-adjoint 
form 


2D (— 0 geil) a 


with functions g,€%,, the formulas for the ¢,’s and the sum )) “/C (e, u)'—*+ 


l=n—j 
simplify somewhat. 


III. Approximation to the solution of the partial differential equation 
We will consider the boundary-initial-value problems (I, 3, 4,5) under the 


additional assumption that the solution (x, ?)€%,. This assumption is con- 
venient for the error estimate we shall give. We use the notations 


ii,= [T flemsys Xif= [X fleas, 


for any {€%, and any point x; in the interval [0,a]. Then we assume 
A1: There is a sequence {N} of integers N, r=1, 2, ..., O< N<N4,, such 
that for all partitions of x,;=ih, h=a/N,, and for all N=N+1—k—1 points 
x; of a partition, =k, k+1,..., N—41, remainder operators of the form 
0, or 
< ” 3 hi t,x C,, i » OF 
Pe} - 25 Oty Xitulion + 2 Biitahien + 2, ™ 
> hwy; Cf 
x=1+s 
exist for any /€%,, where k20 and /20 are fixed integers independent of N.. 
The integers o,, T,, 02, T, may depend on #, the sums o,+ 1, and o,+T, are 
bounded by a fixed integer, and boundary terms h*y;,C,f occur only in 
operators Pf assigned to points x; in a certain neighborhood of the boundary 
points x=0 or x=a. Moreover, the coefficients «; ;4, Bi,i+0) Vi. are € R, and 
all a, ;(#)=-0 on [0, 5]. 
We set a; ;,,,=0 for u>t, and ~<—o, when 1, and 9g, are the values cor- 
responding to 7, and consider the matrix 


Xe ke Xe k+t c++ Oe nNp—1 
a ap k+l *+ 
A= “wa +1,k+ 
2n,—1,k o.6 hoe) O% 0.0 04 Oy,—1,Ne—! 


Then we require 
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A 2: det A+0 on (0, 5]. 
For-the exact solution (x, ¢) the following system of equations holds 
T Tt; 
ps 05 tu lity Mit pe — nh 2 Biitn K+, 


bo, 
, or 
t, y hx 
iN a Say, «xx, OF 
D %,itulitn — WOH + eat Ms, 
B=, n 
po Wn 97, Cys 
x=s+1 
where we have used the relations 
Xian Mitn = — Tipp tisu they 


and the abbreviation 

P£u=h*r(t), 7,ER. 
With the undetermined quantities v,, i=k,...,N—Jl, and with v,=4; for 
t=0,1,...,k—1 and i= N—/+41,..., N, 


pag tit Tin %4+n— Wh 2 Bits Yitn = 
t, | 0, or (IIT, 1) 
oe ee ew n 
preg tite hit , vr p hin-* Vi, x&x» OF p hin yy, 2 xs 
=e] x=s+1 


is a system of N ordinary linear differential equations of order m with a solution 
v;=%(x;, t) on [0, 6] which satisfies the initial conditions 


dh 
aw #(%j,0+) = ui) (%)), &=0,1,....m—1. 


We assume that for i+y<k and for 1+y>N—!/ the functions u(x;,,,,¢) can 
be approximated by the known functions #;,,(¢) which are m-times continuously 
differentiable on [0, 6]. Then we consider the approximate system 


D Mtn Tite ity — wD Biitn Uj4n= 


u=—O; 
. 0, or (III, 2) 
preg tistalisn * } hn" y, Cy, OF >» ia~* Vi, x n> 
x=1 x=s+1 . 


with the undetermined quantities U,, i=k,...,N —l, and with U;=u; for 
t=0,1,...,R—1 and i=N—/+41,...,N. In this system (III, 2) all coeffi- 
cients are known functions €R. For det A +0 this system has a unique solution 
U;=U; (t) satisfying the initial conditions 


dh 
awe Ui(O+) = 4h) (%;), &=0,1,...,.m—1. (III, 3) 


We investigate the differences z;=«u(x;,?)—U;(t) of the solutions of systems 
(III, 4) and (III, 2). They are the solutions of another system of ordinary 
differential equations 


Tt; 


D Mitty liteteu—*” LD Bi itntien = — AO, (III, 4) 


de S| a=-—0O, 
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with z;=(x,,t)—%,(t) for i<k and for i> N —/ and these z,’s € R. The initial 
conditions for (III, 4) are 


a’ 
ar %s(0+) =0, k=0,1,...,.m—1. 


By analogy to (III, 1) being assigned to the given problem, (III, 4) is assigned 
to the homogeneous problem 


aia s=0 
pir 2(%, 0+) =0, 
C,z(x,t) =0 


For zero initial values the system (III, 4) has a uniquely determined solution 
which is still dependent on the unknown function u(x, ¢). If we can estimate 
|z;| we have a measure of the accuracy of our approximation. For the following 
particular class we give such an estimate and sufficient conditions for the con- 
vergence of our method. 


We define the matrices 


Br,r Br,r+1 Hboveser 
B=h-" Bert r Br+i,e+t supe sees 


ore eee ee eee eee By, —1, Nr-> 


F =. oe, k+p Thy (tp +n — Ma+p) 


Terre Te Tree rreyry tr. 
T% 


Re ONn,—1,Ne—1+u Lnp—t+p (MNe—t+y — tia 








= | 
D> Br, rta (“nin — Meta) 
n=—O, 


b=h-"|0 , 


CP atta (Mn Uy,—1+2) 








Zp Tr 


2n,-1 N,-1 


where the number of nonzero components of @ is m, and of D is n,, , and n, 
independent of N. £;;,,=0 for m>t, and for ~<—o,. Then we assume 

A343: The operator T is independent of x. For any of the partitions of A 1 
the matrix A in the corresponding system (III, 4) is independent of ¢. The 


matrix B is given by y(t)B where z(t)€R, and where B is independent of ¢. 
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A and B are symmetric and commutative and all eigenvalues of A4A—B are 
distinct. 
In this case (III, 4) yields 


A(Tz)—yBz=—a—b—hih?"r. (III, 5) 
To solve the corresponding homogeneous system 
A(TZ)—y7BZ=0 (III, 6) 
we set 
h, 


Z=g(t)h, h=| "+ 


’ 


Ay,-1 
where g(t) is a function of ¢, and h is independent of ¢. Separating the variables 
we have a system of algebraic equations for h 
(AA— B)h=0. 
As a consequence of A 3 
det (A A — B) =0 


has real and distinct roots A,, v=1,2,...,.N. Then for each » the equations 


Tg —A,xg=0 (III, 7) 
and the adjoint equations 
T*s—A,xs=0 
have each m linearly independent solutions g,,(t) and s'?—")(t); 1=4, 2,..., m; 
& 3, "ER. Moreover, the systems 
(Ai,— B)h=0 
have N solutions h 
kr 
h hyst» 
hy, 1,» 
which we can assume to be orthonormal. Therefore by SCHWARz’s inequality 
N,-l al 
D> |A,,+| SVN. (III, 8) 
r=k 
Then let the solutions g,; and s!~" be such that 
l Boa Bye oe Som 7 
Y, = _ 7” a alg - 
io le Som 
is a fundamental matrix for the first order system equivalent to (III, 7), and that 
Sy Sy2 s+ Sym | 
-— —- = ... & 
aig) signs)... aber 
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is a fundamental matrix for the adjoint system of this first order system with 
the property 


a 4: 0),. 2a1@ 
bx sii} ten ‘stg! Som 01 : 
* = 
Yy Py = . = " *. eZ 
-1 =1 ° ° 
"5. sii not Lo. gl gi 1 0 
7=0 j7=0 seentauns 0 = 0 1 


where s=s,,;, g=g,,. It is always possible to find fundamental matrices 9, 
and y, [8]. 
For equations (III, 6) the vectors 


Rhg,;, v=1,2,....N, %#=4,2,....,% 
are Nm linearly independent solutions: for if there exist constants k,; (not all 
zero) such that 
Dh My Bri = 0, 
»,4 


then 
Lh, g, =0 
with not all =) k,; £,4==0 because not all k,; are zero and the g,,;’s for a fixed 


+ 
vy are linearly independent. This is a contradiction; for the h,’s are linearly 


independent. 
Therefore 
h, byi 
@ _ ” by 
h, qr 


is a fundamental matrix for the first order system equivalent to (III, 6) and 
similarly 

h, s,; 

h it 


rae 


h sim—1) 


vt 


is a fundamental matrix for the corresponding adjoint system. We find 


78 ..'G 
m—1 0 1 
v+o—|S stontneft|=|. 0. 2]. 
r=0 1 0 
a 


Thus the solution 2 of (III, 5) with zero initial values can be written 


= — [h, g,;(t)] f [sip (x) hE (a(x) + B(x) + he-* r(x) ] dt 


or 


2=— Dh,z,:(t ) fsfa-(2 ) (h¥ a(t) + h¥ b(t) + he" h¥ r(t)) dr 
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We define 
similarly for ||a||, nh oak 
ll =,_,mag (1D 


oak 2,...,N 


Then, because of ||h,|] <1 and (III, 8) 
[el] SPI] (0 |||] + |] 8|] + 4°-* YN || r|I) >| fees si '(r)dz}. (III, 9} 


The functions g,,(¢) and s!~"(#) depend on 4,, are continuous on [0, 6], and 
are therefore bounded for each particular A,. By P 3 


jim |7,(4)| =0, 
further 
YNhs ja. 
We assume 
A 4: For the partitions considered in A 1 


[A] |e], =0, jim ||h,|]|]5]| =0, Jim re-*~# ||, |] | |] = 


h—0 jim | h—>0 


and 
>I fet slvr) dt .f, 


where J is a uniform bound. 
Then on account of (ITI, 9) 


lel] S[lAoll ellal| + me|| bl] +4-*-F Yall) 2 (LIL, 10) 


and our procedure is uniformly convergent for the partitions of A 1. 


IV. Applications 


We now consider some simple but important examples to illustrate our 
methods and the principle of the convergence test. 


1. Given the equation 


(7) 
a ST +H(xt), 


OsxS1i, OStsSb, c*>O0, c*constant, /ER, 


with the initial condition 
u(x,0)=u*(x), O<%*<1, 
the boundary conditions ve 
4 (0, t) = u(t), 
u(1,t)=4, (4), OStsSd, 


which we assume has a unique solution «(x,?)€%,. For partitions %,=th, 
h=1/N,, N>1, i1=0,1,..., N, a set of remainder operators P? like those re- 
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quired in A 14 is given by 


P? f =h* fi! — (fa — 2h + har), 4+=1,2,...,N, 
with 


h* 
PP] S az llAll- 


In this case we have kR=/=1, N= N,—1, o,=1,=0, 0g=1,.=1, constant coef- 
ficients &; :44, B:i+., and «,,;=1. Det A=1 and thus A 2 is satisfied. System 
(III, 2) is given by 
dU; : 
— = (Ga— 2U,+ Us) =f, 
with Uj=%), Uy.;=%, and initial conditions U;(0)=«*(x,). System (III, 4) 
becomes 


dz; c? 
a; — >} (4-1 — 24,+2,41) =he'x,, 





with z)=%y,,=0. The matrix B=B is symmetric, A and B are commutative, 
the operator T is independent of x, the eigenvalues of (AA—B) are given by 


i= 


and are all distinct. Thus A 3 is satisfied. We find 


2c? 
i (1—cosp,), p~=2ahy, v=1,2,...,N, 





81 = &, = exp[— At], 
S,1 = S, = exp [A, t] , 
h, , = V2hsin (yu, i), 


||”, || < 2h, 





and 
h 
rl] < * |u]). 
Further, 
y fel s(t) dz —_ ¥ 2 2 (4 — g,(t)’ 
wals S26 (1 — cos p,) 





1 1 
< > : r — ‘ —s 
=] == a= =—_ ee 
’ a(t a) a 6c (1 a) 





Thus A 4 is satisfied, from (III, 10) it follows that 


h2 2 [|u| 
l|] see 


and our approximate solutions converge uniformly as h? toward the exact solution. 
2. We next consider the same problem as in example 1 but with boundary 


conditions 
4 (0, t) = u(t), 


a _ 
oy H(t. 4) =m (2), Oositsd. 





| 
| 








a 


en ett 


we me 
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We-use the remainder operators P? and the relation 


, hh? 
— t+ twa = — twa — fn 
valid for any {€8,, where /y,, denotes the second derivative of { with respect 
to x at a point x, %y<%<4%y4,. Thus we write 


Lad 


7 —_? 2 
Py 4 = h® uy — (Uy—1 — My) + hi + = 4N+1 


to match the boundary condition at x=1. System (III, 2) is given by 


dU; _ 
dt 





2 
a (Ua — 2U,+ Us) =f; 
with 
Uy = ty,— Uy + Una = hi, 
and with initial conditions U;(0)=«*(x;). System (III, 4) becomes 


az; c® 
at << (2-1 — 22, + 241) =he*y,, 





with 20> 0, = 2y t+ 2y41 => uUN4+a ° 
The constant matrices A and B are again symmetric and commutative, and 
the eigenvalues of (AA — B) are 
2v—1 


A,= 73 (1 —cosy,), My = NNT? y= 4,2,...,N, 





and they are all distinct. Thus A 3 is satisfied. We find 


81 = & = exp[— A, ¢), 
S,1 = Ss, = exp[A,Z], 


h,, = 2y% sin (u, 2) 
2+(- 1)” hk sin [u,(N + 1)] cos & 








with 
lim || h,|| =0 
h—0 


for ail y values. Further 











N N he ‘ 
De J sl s, (7) dt < 2 2c*(1 — COS py) < 6c? (1- 7 ’ 
2 


and 
2 
llal|=0, [J <sS][u"]], m=4. 


Thus A 4 is satisfied. From (ITI, 10) it follows that 


0 al 1 
Nell S11 (Il + lll) ssa 





and our procedure is uniformly convergent. 
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3. We again consider the same problem as in example 1 but now we use 
higher order approximations. Let o=5, u(x,t)€B,, and let the remainder 
operators P;> be given by 

PP f= 1 (fia + 10f;' + fia) —12(fKi— 2h +h), 1=1,2,...,N 
with 
[Pfs * ||]. 


We have k=/=1, o,=11=0,=1,=1, and constant coefficients «; j4,, Bj j+n- 
Using the recursion formula 


det Ay = 10 det Ay_; — det Ay-_»; 


it can be shown that det4210. Thus A 2 is satisfied. System (III, 2) then 
becomes 





Hi ima 4 49. a0 4. Bits _ 120 (U1 — 2U, + Usa) = fa + 10f, + hiss 


with 
ain a ——_ wf 
= aie Wart = ah 








and with initial conditions U,(0)=«*(x,). System (III, 4) becomes 
a 








=i —1 + 140 + ts a? = 2 Z-41+ 102; +z 2541) = h'c*r,, 





dz d 
fe Mya, ae ==(Q. 


Both matrices A and B are symmetric and commutative, and the eigenvalues 
of (AA—B) are 














_ 12c? 1—coshay mm 
A, = h? 5+coshay’ y=1,2,...,N, 
and they all are distinct. Thus A 3 is satisfied. We find 
&, = €xp i- A, t] , S, = exp [A, t] , 
h,, = \2hsin(xhyvi), 
and 
N,¢ N 
J h?(5+ cosh xv) — es TE 
Da J sl s,(t) dt ard Da 12c?(1— cosh xv) (1 é )< 6c2(1- =) ad 
12 
Thus A 4 is satisfied and 


| «6% |] . 


ll-ll s mee say 
4. Given the equation 
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with the initial conditions 
u(x,0) = u*(x), £u(x, 0+) =ui)(x), O<%*<a, 
the boundary conditions 


(0, t) =t(t), 
u(a,t)=a,(t), OStSd, 


which we assume has a unique solution u(x, t)€ 8,. We use the same remainder 
operators as in example 1. System (III, 2) is given by 


BU; 1 , 
an — <5 (Ga — 2U;+ U, 41) =0, $m 1,2,...,N, 
system (III, 4) is given by 


d*z; 1 
aa ee h? (2, i-1 22; +2 +41) = vit 


29 = 2y41= 0, 


and the eigenvalues of (4A — B) are again 


== (1—cosahv), v=1,2,...,N. 








We find 

81 = Ars cos Vt, 8,3= A>! sin VA, 
sl) = — Arisin V/A, ¢, patna 
h,, = 2hsin(xhyvi), IIrI|<- lH), 

and 

2 >| fen (t) s@(zr)dt} < ye [1 — cos YR] $29 45+ < _ +. 
Thus 
Wey2_ yay 

llz||s = 3(12— mm) -|| ||. 
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Handbook for automatic computation 


Preparation of a handbook for automatic computation, in five or more volumes, is now 
under way for publication by Springer-Verlag. It will appear in the series, ,,Grundlehren 
der Mathematischen Wissenschaften“ (‘Yellow Series‘). Editors are 


F. L. BAvER, Mainz 

A. S. HOUSEHOLDER, Oak Ridge 
F. W. J. OLvER, Teddington 

H. RUTISHAUSER, Zurich 

K. SAMELSON, Mainz 

R. SAUER, MUNICH 

E. STIEFEL, ZURICH 


The purpose of the handbook is to provide a collection of tested algorithms for mathematical 
computations of all sorts: the solution of finite and of functional equations, methods of 
approximating functions, the evaluation of special functions, etc. These algorithms are to 
be written in ALGOL, hence will be usable on any machine for which a suitable translator 
is available, and even without a translator can be used as a model for programming. It is 
evident that such a collection could have no general utility unless written in some common 
program language. The descriptive language will be English. 


As plans now stand, the organization of the series will be as follows: Volume 1a will contain 
a description of the use of ALGOL, and Volume 1b a description of the structure of trans- 
lators. These introductory volumes are the only ones that will not be made up primarily of 
actual algorithms. Volume 2 will be devoted to the solution of finite equations, linear and 
nonlinear, including the determination of characteristic values and vectors of matrices. 
Volume 3 will be on functional equations, especially differential equations, ordinary and 
partial, and integral equations. Volume 4 is concerned with méthods of approximation, and 
Volume § the evaluation of particular functions. It is possible that certain algorithms, such 
as those for solving inequalities, for mathematical programming, for statistical computations, 
and the like, that do not seem to fall naturally in any of these areas, may be reserved for 
a sixth volume. Each algorithm is to be accompanied by enough explanatory information 
to make it understandable, along with whatever information is available on speed, accuracy, 
range, or, more generally, for judging the effectiveness of the algorithm for a given type of 
problem. In any event, only pretested algorithms will be published. 


Before the appearance of the volumes themselves, the algorithms will be prepublished in a 
series of supplements to the journal, ,,Numerische Mathematik". This is partly to make 
generally available each algorithm at the earliest possible time. But in addition to this, it 
provides the possibility for including in the handbook itself additional information, and even 
corrections, th t might come in from users. 


Contributions are earnestly solicited. For the present, at least, these must necessarily be in 
the form of actual algorithms, along with information as to the extent and mode of testing 
the algorithm, estimates of accuracy, and experience in using it. Untested algorithms will 
not necessarily be rejected ipso facto, but their use must necessarily await actual test. As 
algorithms are published, information relating to published algorithms also will be welcomed. 
Contributions may be sent to any of the editors named above. 
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